ncar/topo

NCAR Global Model Topography Generation Software for Unstructured Grids

topo's Issues

regridding FV3 c540/CONUS refinement factor gives isolated zeros

When a refinement factor on the c540/CONUS SCRIP file is regridded to the nc3000 grid in cube_to_target.F90, isolated stripes of ZERO are produced. The field on the SCRIP is fine, i.e., minimum value=1.0.

-z option not working

For some reason the command below doesn't result in ldevelopment_diags=.TRUE.

./cube_to_target --grid_descriptor_file='../../regression-test-data/ne30pg3.nc' --intermediate_cs_name='../../regression-test-data/gmted2010_bedmachine-ncube0540.nc' --output_grid=$ogrid --coarse_radius=$Co --fine_radius=$Fi -r -u '[email protected]' -q 'output/' -z

The additional output isn't written.

If I put the '-z' option right after the '-r' option I get the error:

"Required argument not specified: name_email_of_creator"

Add script for energy spectra

Create DOI with Zenodo once we are ready with new release

https://zenodo.org/account/settings/github/

add a file name modifier when running w/ no leak (or vice versa)

topo file names should have some indication whether no leak is on/off

Mars

Source topo dataset for Mars

/glade/p/cgd/amp/pel/topo/rawdata/mars-rawdata.nc

is missing data North/South of 88.2N/S which trips error checks in bin_to_cube.

ncdump below:

lat = -88.1949168841954, -88.1871353846802, -88.179353885165, ...

ESMF grid descriptor files for spectral-element refined grids do not have refinement factor variable

/glade/p/cesmdata/inputdata/share/meshes/ne0ARCTICGRISne30x8_ESMFmesh_c20200730.nc

does not have refinement factor variable!

Note: in ESMF format we have unofficially decided that the refinement variable is named elementRefinementRatio
(this is what is used on the MPAS regional refinement grid descriptor files)

Note that the SCRIP version in the CAM coords directory does have refinement factor variable (rrfac):

/glade/p/cesmdata/inputdata/atm/cam/coords/ne0ARCTICGRISne30x8_scrip_c191209.nc

add fast (low res?) regression tests

A couple of minor changes

1. remove (Nsw, Nrs) "_Fi001" from file name.
2. make --fine_radius optional (default 0)
3. introduce develop option that turns on writing intermediate ".dat" files, writing ridge variables that are not used by CAM (currently they are on the topo file - remove), print statements and other files used for development and not production.

Adjust Smoothing to Match MPAS smoothing

An issue to discuss updating the software to match the smoothing of MPAS.

implement more accurate estimates of jmax_segments when using variable resolution grids

dealloc error for non var-res grids

I just ran topowopts with this command for a ne30pg3 grid:

./cube_to_target --grid_descriptor_file='/glade/work/aherring/grids/uniform-res/ne30np4.pg3/grids/ne30pg3_scrip_170611.nc' --intermediate_cs_name='/glade/p/cgd/amp/pel/topo/cubedata/gmted2010_modis_bedmahcine-ncube3000-220518.nc' --output_grid='ne30pg3' -r -d -m --smoothing_scale=100.0 -u "Adam Herrington, [email protected]"

After almost completing this job (it did get past writing the final file), I got this error:

 writing GBXAR  data   2.1727097902297268E-004   3.0458649516175928E-004
 done writing GBXAR data
 close file
At line 973 of file cube_to_target.F90
Fortran runtime error: Attempt to DEALLOCATE unallocated 'target_rrfac'

I suspect the target_rrfac pointer is not being allocated in read_target_grid when there is no rrfac on the SCRIP file.

Do we need --zero_negative_peaks and --ridge2tiles namelist options?

just checking ...

lstop_after_smoothing=.true. and refined grids

lstop_after_smoothing=.true. will not work for refined grids since rrfac is not read in ...

Re-enable regional refinement for topowopts branch and other ridge finding mods

I've merged (on my fork) my latest mods with the up-to-date NCAR/Topo/topowopts. There are some issues that we should discuss:

1. I'd like to eliminate use_prefilter in favor of just keying off the fine_radius parameter. Code is now logically consistent with that.

2. I'd like to eliminate separate calls to several subroutines that depend on the value of lregional_refinement flag. I believe we should just let the codes work with rrfac=1 everywhere if refinement is not desired. The case lregional_refinement=.false. should just trigger the creation of an rrfac array that is =1 everywhere.

Issue with not smoothing over LANDFRAC=0 areas on high-resolution intermediate cubed-sphere grid

Plot below shows the difference between smoothing everywhere and not smoothing when LANDFRAC=0 with ncube=3000 intermediate cubed-sphere grid (~3km):

and ncube=540 (~12km)

There are significant differences inland. Why?

The problem: Binning from ~1km LANDFRAC to 3km cubed-sphere picks up on isolated "inland" ocean points that are not picked up when using ncube=540 intermediate cubed-sphere grid. See this regional plots over Himalayas of LANDFRAC when using ncube=3000:

and ncube=540:

There are sporadic LANDFRAC=0 points (dark blue points) where the Laplacian smoother will not be applied in the ncube=3000 LANDFRAC but close to none in ncube=540 LANDFRAC.

ESMF grid descriptor format does not work for lat-lon grids

There is no way to tell that a grid is lat-lon from ESMF grid descriptor format. Probably need option in code to tell it is lat-lon so that the final topo file is 2D.

Sample ESMF lat-lon file located here:

/glade/p/cesmdata/inputdata/share/meshes/fv0.9x1.25_141008_ESMFmesh.nc

Writing intermediate smoothed topo file not working for Laplacian smoother

errors out for uniform resolution SCRIP files that contain rrfac

The code currently turns on regional refinement if there is an rrfac in the SCRIP file, but some uniform resolution SCRIP files may contain rrfac, causing the software to fail for those grids.

Propose setting rrfac_max = 1 for such a uniform resolution grid, to cause the program to ignore rrfac in the SCRIP file and otherwise carry on as if it were a uniform resolution grid. The code currently will error out if you set this to 1 and rrfac is in the SCRIP file.

Here's a rotated ne30 SCRIP file that contains rrfac, for testing:

/glade/scratch/patc/AMP_REPO/ne0np4.SAMwrf01.ne30x1/grids/SAMwrf01_ne30x1_np4_SCRIP.nc

@patcal

re-introduce LANDFRAC

lzero_out_ocean_point_phis=.true. algorithm needs LANDFRAC

Re-introduce LANDFRAC on intermediate cubed-sphere file, in cubed_to_target.F90 but don't write LANDFRAC to final topo file.

Output file name character issue

From standard output:

Final output file :: ./output/°uìÍxHîjiî@ üp üjiîü°
(jiî<<Hî 
c;HîjiîPj?Hî``xz2ìjiîj?Hîxz2ìjiîx01ì(  üàmiîjiî\JHîüji 
ü 
ü 
ümiîP1k÷ë1ìjiî``ýãì`zsJìXt

-v, --rrfac_manipulation

@adamrher suggested that "-v" should just turn on smoothing of rrfac as well as integer operations on rrfac (needed to spectral-element grids) but that even for MPAS it might be good to bound rrfac between 1 and rrfac_max.

Add regression test to ensure bit-for-bit changes before making changes

Before we make any major changes, creating a regression test to ensure we get the same results will be beneficial.

This issue serves as an area to discuss creating this regression test.

Always read in PHIS on GLL grid for spectral-elements

Decide on naming convention for PHIS on GLL grid when using physics grid option. This will require changes in dyn_comp.F90 on cam_development (i.e. coordination needed).

See also issue ESCOMP/CAM#552

Scientifically validate high uniform and high resolution variable resolution MPAS topo files

I have been generating several MPAS topo files that need to be scientifically validated:

1. currently doing rrfac manipulation. Should we? (subtract uniform resolution topo files with resolution of low res and high res areas; sanity check)
2. when MPAS reaches ~3km in refinement region; is the topo "well-behaved"? Would need to do higher resolution intermediate cubed-sphere resolution and check sensitivity to that)

Command used for 60-3km grid:

./cube_to_target --grid_descriptor_file='/glade/work/skamaroc/MPAS_meshes/CESM_MPAS/60-3km/West_Pacific/mpas_x20.835586.wpac.1_esmf.nc' --intermediate_cs_name='/glade/p/cgd/amp/pel/topo/cubedata/gmted2010_modis_bedmahcine-ncube3000-220518.nc' --output_grid='mpas_x20.835586.wpac' -u 'Peter Hjort Lauritzen, pel@ucar' -q 'output/' --smoothing_scale=50.0 -r -y 20 -v -m

-m = smoothing over ocean points
-v = rrfac manipulation
-y = refinement factor

Command used for 60-10km grid:

./cube_to_target --grid_descriptor_file='/glade/work/skamaroc/MPAS_meshes/CESM_MPAS/60-10km/mpas_x6.999426.wpac.1_esmf.nc' --intermediate_cs_name='/glade/p/cgd/amp/pel/topo/cubedata/gmted2010_modis_bedmahcine-ncube3000-220518.nc' --output_grid='mpas_x6.999426.wpac' -u 'Peter Hjort Lauritzen, pel@ucar' -q 'output/' --smoothing_scale=50.0 -r -y 6 -v -m

Topo files here:

/glade/p/cgd/amp/pel/topo/grids/mpas-60-3km-WestPacific/mpas_x20.835586.wpac_gmted2010_bedmachine_nc3000_Laplace0050_20220520.nc
/glade/p/cgd/amp/pel/topo/grids/mpas-60-10km-WestPacific/mpas_x6.999426.wpac_gmted2010_bedmachine_nc3000_Laplace0050_20220520.nc

merge TopoCESM2 to master

Hi @gold2718

@PeterHjortLauritzen and I talked about merging the TopoCESM2 branch to the head of master. There have been many updates to TopoCESM2 so I think it's time to just bring it in to master.

I'm making some tweaks right now and so will commit a final version of TopoCESM2 in a couple hours, but I just wanted to put out a feeler.

bug fix for variable resolution configurations

I noticed a bug in my implementation of the regional refinement enhancements in branch TopoCESM2. In particular, the refinement factor should have the limiter applied to it before being called by the topography smoother. The refinement factor is instead going through the limiter step after the topography smoother is called. This means that the smoother isn't calling the correct refinement factor, which is a bug.

This bug has been fixed in a forked-off branch TopoCESM2.2 in my git repo:

https://github.com/adamrher/Topo/tree/TopoCESM2.2

I also updated the experiment_settings.make file by grouping the three (3) supported variable resolution grids (ARCTIC, ARCTICGRIS & CONUS) together under the sub-header "functionally supported vr-grids." I also had switch the ARCTIC & CONUS settings to default to including anisotropic gravity wave fields, as req'd by CAM6 physics.

Updating to bedmachine landfrac

Rene Wijngaard has updated our MODIS land fraction dataset to GrIS and AIS high resolution bedmachine data over Greenland and Antartica. The difference between updated dataset and GMTED2010 is here:

Interpolation to grids with <= 3 km grid spacing produces too smooth of height field

When interpolating to grids that with 3 km grid spacing, the topo tool produces a height field which is significantly lower resolution than what is expected. I belive this is caused by the double interpolation from the ~1 km GMTED2010 data to the 3 km cube sphere then to the target grid. Then, with ridge analysis the height field is smoothed even further.

Below are two plots of the height field over conus. The first, is the x5.6488066.grid.nc (15-3 km variable resolution mesh) with the 3 km refinement centered over CONUS with the height field that was created with this tool. (I believe the height field outside of the continent lines is due to a lack of landmask).

The second is GMTED2010 dataset interpolated to the x1.65536002.grid.nc (3 km quasi-uniform grid) gird using the MPAS-Model`s init_atmosphere core.

These files can be found on glade at: /glade/scratch/mcurry/cam-mpas-topo-compare

Is there a way interpolate GMTED2010 to a higher resolution cube sphere?

lfind_ridges=.true. does not work with Laplacian smoother

Validating new topography dataset and smoother

[transferring text from issue https://github.com/NCAR/amwg_dev/issues/116 to this repo as it is relevant; the plots in this post are using GMTED2010 elevation data only (no bedmachine)]

New topography file

We have produced a new topography file where we use a Laplacian smoother (instead of the default distance weighted smoother) and turn off smoothing where the land fraction is 0 (i.e. no smoothing over ocean points; partial ocean points are still smoothed). We are curious if not smoothing over ocean points has an impact in coupled model setup since the temperatures passed to the ocean, where previously the topography was "leaking into the ocean", could be biased cold (assuming temperature decreasing with height on average).

The Figure below shows energy spectra for different topographies:

A couple of observations regarding the energy spectra:

• The energy spectra is the same when using a 3km intermediate cubed-sphere field (ncube3000, orange line) compared to using a 16km intermediate cubed-sphere grid (ncube540, red line) when computing energy spectra on 1 degree finite-volume grid. That means, we can assess smoothing using the ncube540 intermediate cubed-sphere grid which is much more computationally efficient to run the smoother on compared to the ncube3000 intermediate cubed-sphere grid.

• Below show difference between smoothing everywhere but using the ncube=3000 grid versus ncube =540 grid:

Basically this plot verifies that smoothing on ncube=3000 and ncube=540 using the same physical smoothing scale leads to almost the same result (as one would expect). For comparison distance weighted smoother minus Laplacian smoother (PHIS) below:

[very few differences!]

• Interpolation matters for computing energy spectra: bi-linear interpolation from 1 degree spectral-element grid to 1 degree finite-volume grid produces damping (compare blue and red lines). Please note that the two grid "see" the exact same smoothed topography on the intermediate cubed-sphere grid. Hence when comparing different levels of smoothing always use the same target grid.
• The energy spectra with new and old smoother are very similar (compare red and green lines).
• Not smoothing over ocean does add some energy at many scales (compared cyan and red).
• Black line is unsmoothed topography (for reference)

The next two plots show Greenland area and South America topographies (all mapped to fv 1 degree grid):

• ncube540 raw: raw topography on ~16km intermediate cubed-sphere (ncube540)
• CESM2: Current default topography
• nc540: Laplacian smoothed topography (nu=28E7) - old CESM1.5 smoothing level when the SE topography was smoothed using SE (and not in the topography software as done now; which is needed for the ridge scheme variables)
• nc540;Co11: Laplacian smoothed topography (nu=20E7) that matches C011 smoothing (default CESM2) on the ncube540 grid (for the ncube3000 grid that is equivalent to C060 smoothing radius)
• nc540; no-ocn: Laplacian smoothed topography (nu=20E7) that matches C011 smoothing (default CESM2) but no smoothing over ocean.

The cross section plots (lower right) are though 30S and 65N.

The following observations are made:

• Smoothing over peaks seems to be roughly the same for all smoothing options -> not smoothing over ocean does not seem to affect peaks a couple of grid-points inland.
• One can visually see no smoothing over ocean.

Some mathematical details on the smoothing operator:

Support for ESMF (NUPC,CMPS) grid descriptor file format

The SCRIP files use the following naming convections (example from MPAS):

netcdf mpas_scrip_60-3_andes {
dimensions:
grid_size = 835586 ;
grid_corners = 7 ;
grid_rank = 1 ;
variables:
double grid_area(grid_size) ;
grid_area:units = "radians^2" ;
grid_area:long_name = "area weights" ;
double rrfac(grid_size) ;
rrfac:units = "dimensionless" ;
rrfac:long_name = "normalized dc" ;
double grid_center_lat(grid_size) ;
grid_center_lat:units = "radians" ;
double grid_center_lon(grid_size) ;
grid_center_lon:units = "radians" ;
double grid_corner_lon(grid_size, grid_corners) ;
grid_corner_lon:units = "radians" ;
double grid_corner_lat(grid_size, grid_corners) ;
grid_corner_lat:units = "radians" ;
int grid_imask(grid_size) ;
int grid_dims(grid_rank) ;

The ESMF version is:

netcdf mpas_esmf_60-3_andes {
dimensions:
nodeCount = 1671168 ;
elementCount = 835586 ;
maxNodePElement = 7 ;
coordDim = 2 ;
variables:
double nodeCoords(nodeCount, coordDim) ;
nodeCoords:units = "degrees" ;
int elementConn(elementCount, maxNodePElement) ;
elementConn:long_name = "Node Indices that define the element connectivity" ;
elementConn:_FillValue = -1 ;
byte numElementConn(elementCount) ;
numElementConn:long_name = "Number of nodes per element" ;
double centerCoords(elementCount, coordDim) ;
centerCoords:units = "degrees" ;
double elementArea(elementCount) ;
elementArea:units = "radians^2" ;
elementArea:long_name = "area weights" ;
double elementRefinementRatio(elementCount) ;
elementRefinementRatio:units = "dimensionless" ;
elementRefinementRatio:long_name = "normalized dc" ;
int elementMask(elementCount) ;

jmax_segments

For low resolution grids (e.g. ne16, mpas480) source code change is needed:

replace

jmax_segments = MIN( jmax_segments, 10000 )!for low resolution grids this line needs to be uncommented (makes code slow!)

with

jmax_segments = MIN( jmax_segments,100000 )!for low resolution grids, e.g., mpasa480

in cube_to_target.F90.

Note: if jmax_segments is large then high resolution grids make the model run very slow ...!

Maybe jmax_segments could be optional argument?

new 30arcsec GMTED dataset

New topography dataset. It is the old GMTED2010 dataset merged with the GrIS and AIS high resolution bedmachine data, for more accurate representation of ice sheets. Here is what a ne30pg3 topofile looks like using this new dataset minus the old dataset:

Among other changes, the ice sheets have steeper margins, which seems to produce larger subgrid heights. This is interesting b/c if you recall, SGH30 was artificially increased over Greenland in the old dataset to bring the surface winds into better agreement with observations. @whlipscomb @lofverstrom

Once #9 is merged, I can modify the Makefile to point to the new 3km intermediate cubed-sphere dataset.

failure to read command line arguments with Intel compiler

Remove LANDFRAC from topo software

Remove landfrac from topo software as it is not read in by CAM. It is a different landfrac than used in the model, so it's really just a point of confusion and should be removed.

can we remove landm_coslat? or at least optionally remove it from topo file?

Add GLL grid topo as well as lat_d, lon_d to topography file.

I think we should have topography on both the phsygrid and the GLL grid in topography files for physgrid configurations e.g., ne30pg3

refinement factor (rr_fac): introduce option to turn off/on internal manipulation of rr_fac

Code for manipulating rr_fac is tailored for SE.
MPAS does not need extra smoothing of rr_fac.

Hence make manipulations optional.

argument order issue

./cube_to_target --smoothing_over_ocean --grid_descriptor_file='/glade/p/cesmdata/cseg/inputdata/share/meshes/ne16np4_ESMFmesh_cdf5_c20211018.nc' --intermediate_cs_name='/glade/p/cgd/amp/pel/topo/cubedata/mars-ncube3000.nc' --output_grid='mars_ne16np4' --smoothing_scale=200.0 -u 'Peter Hjort Lauritzen, [email protected]' -q 'output/'

works but when moving --smoothing_over_ocean to the end

./cube_to_target --grid_descriptor_file='/glade/p/cesmdata/cseg/inputdata/share/meshes/ne16np4_ESMFmesh_cdf5_c20211018.nc' --intermediate_cs_name='/glade/p/cgd/amp/pel/topo/cubedata/mars-ncube3000.nc' --output_grid='mars_ne16np4' --smoothing_scale=200.0 -u 'Peter Hjort Lauritzen, [email protected]' -q 'output/' --smoothing_over_ocean

the model just hangs!

Consistent representation of sub-grid ridges for variable resolution grids

@JulioTBacmeister is working to enhance the topography software to provide consistent sub-grid representation of ridges for variable resolution grids. The ridge information is used by the anisotropic orographic gravity wave drag scheme in CAM6, and it is important for this to be done correctly in variable resolution configurations as to not "double count" owing to both resolved and parameterized orographic gravity waves. This is a high priority task due to the need to provide consistent, correct topography boundary conditions for the roll out of our variable resolution grids in the next cesm release (cesm2.2). I am opening up this thread to discuss progress and/or difficulties in completing this task.

For debugging: add refinement factor to topo file (as an option)