Skip to content

Post process geogrid output to compute GWD input fields directly from source data#281

Open
islas wants to merge 27 commits into
wrf-model:developfrom
islas:compute_gwdo
Open

Post process geogrid output to compute GWD input fields directly from source data#281
islas wants to merge 27 commits into
wrf-model:developfrom
islas:compute_gwdo

Conversation

@islas

@islas islas commented Feb 4, 2026

Copy link
Copy Markdown
Contributor

These changes calculate the necessary fields in a manner that accounts for subgrid-scale orography. This in in support of subgrid orographic parameterization revisions for GWD in WRF.

Currently the gravity wave drag (GWD) input fields from WPS/geogrid are computed from averaging of the source data then performing the calculations for oa1-oa4, ol1-ol4, var, and con. This utility program instead calculates these orographic statistics directly from the source data at a subgrid-scale assuming the domain dx/dy are larger than the source data extent per value.

The code is designed to be run using the entry point ./util/compute_gwdo.py and will use ./namelist.wps if no namelist is provided.

The namelist is read to determine the geogrid domains output, source data tiles locations, and domain spacing. To avoid recalculating lat/lon positioning, the values are read from the mass staggered lat/lon locations within geo_em.d00.nc, geo_em.d01.nc, and so on. Source data is loaded as tiles are needed to compute the orographic statistics for a given domain cell based on the dx/dy extent of the cell at that lat/lon.

Once all calculations are done for all cells within a domain file, that file's original oa1-oa4, ol1-ol4, var, and con are modified with the newly calculated subgrid-scale statistics. Additionally, a MAX_EL field is created. This is repeated for all domains listed within the namelist.

@islas islas changed the base branch from master to develop February 4, 2026 23:35
@mgduda mgduda self-requested a review May 26, 2026 20:30

@mgduda mgduda left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I'll take more time to try out domains over different parts of the globe, but for now I thought I'd leave some minor initial comments.

Comment thread util/compute_gwdo.py Outdated


# https://stackoverflow.com/a/34325723
def progerss_bar( iteration, total, prefix="", suffix="", precision=1, length=100, fill="*", empty="-", end="\r"):

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Would it be better to name this function progress_bar?

Comment thread util/compute_gwdo.py Outdated
Index slices, multiple key-value pairs on the same line,
and other features are NOT supported. For full support
please utilize f90nml.
Spaces between key and value and commas does not matter.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

"Spaces" -> "do not matter"

Comment thread util/OrographyStats.py Outdated
self.box = box

self.mean = np.mean( self.box )
# var (actulally stddev)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

"actulally" -> "actually"

Comment thread util/SourceData.py Outdated

class SourceData:
def __init__( self, name, path ):
self._earth_radius = 6371229.0

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I think the WPS and WRF use 6370000.0 m as the sphere radius (MPAS uses 6371229.0 m).

Comment thread util/SourceData.py Outdated
self._subgrid_m_dx = 2.0 * np.pi * self._earth_radius / self._npts_x

self._tile_data = TileData(
int( int( 360.0 / self._index.dx ) / self._index.tile_x ),

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Rather than repeating int( 360.0 / self._index.dx ), would it be better to use self._npts_x?

Comment thread util/SourceData.py Outdated

self._tile_data = TileData(
int( int( 360.0 / self._index.dx ) / self._index.tile_x ),
int( int( 180.0 / self._index.dy ) / self._index.tile_y ),

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Would it make sense to define self._npts_y as we do for self._npts_x, and to use that value here?

Comment thread util/compute_gwdo.py Outdated

def main():
parser = argparse.ArgumentParser(
"./compute_gwdo.py",

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

What do you think about removing the ./ prefix from the program name? Since this resides in the ./util subdirectory, my guess is that most users will run this as ./util/compute_gwdo.py in practice; but it's possible that users may also have an entirely separate run directory from which they will run this script.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

What about something like <path to script>/compute_gwdo.py?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I think it would be reasonable to assume (expect?) that users have some experience working in a shell, and they will understand that they can either (1) provide a relative or absolute path to the compute_gwdo.py program when running it, or (2) add the path to compute_gwdo.py to their PATH environment variable and subsequently run compute_gwdo.py with no path specified.

Comment thread util/compute_gwdo.py Outdated
)
parser.add_argument(
"-d", "--dataset",
help="Topographic dataset to use, in geogrid format",

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Should we add more info here (or elsewhere?) to indicate that only global datasets on a lat-lon grid will work? For example, using the SRTM 3s data with -d SRTM_topo_3s will fail.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

A quick investigation found that taking a subset will require further work and analysis. For the sake of time I think this is sufficient to at least allow working with the global gmted data and noting it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants