Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

How to reduce topography over a certain region using the read_wrf_nc code

This post was from a previous version of the WRF&MPAS-A Support Forum. New replies have been disabled and if you have follow up questions related to this post, then please start a new thread from the forum home page.

I want to reduce the topography in the geo_em.d02.nc file using the read_wrf_nc over a specified region.

I added the following lines in the USER_CODE.

Code:
 elseif ( var == 'HGT_M') then !
     do xx=1,dim1
      do yy=1,dim2
        data_real(xx,yy,1) = data_real(xx,yy,1)*0.5
      enddo
     enddo

If I do this, the topography of the entire domain is reduced in half.

I want to reduce the topography say, from 120-121.5E and 12-15N only. How should I specify the region in the read_wrf_nc code?
The geo_em file is not in standard grid format. I also do not know how to convert it to lat-lon values.

I'll appreciate any help.
Lyndz
 
I want to reduce the topography (HGT_M) by 1/2 of the geo_em files from geogrid over a specified region.

I added the following line in the USER_CODE of the read_wrf_nc.f

Code:
 elseif ( var == 'HGT_M') then !
     do xx=1,dim1
      do yy=1,dim2
        data_real(xx,yy,1) = data_real(xx,yy,1)*0.5
      enddo
     enddo

However, this reduces the topography over the entire domain by 1/2.
How do I specify the region where I will apply the above code?

I'll appreciate any help on this,
Lyndz
 
Hi Lyndz,
It looks like you're on the right track, and just need to make modifications to specify the exact area for the change. Take a look at this forum topic:
https://forum.mmm.ucar.edu/phpBB3/viewtopic.php?f=45&t=8413&p=14559#p14559
See my response from Sept. 26th, providing an example "if" test for something similar. The difference is that in the other topic, they are modifying the WRF code, as opposed to WPS, and therefore the naming convention for i,j will likely be a bit different. The example shows mods for distributed memory code, hence the "p" (for 'patch') in "ips/jps." Let me know if this helps.
 
Dear Kwerner,

Thank you for this response.

I saw this line:

Code:
IF ( ( ips <= 60 ) .and. ( ipe >= 60 ) .and. &
( jps <= 70 ) .and. ( jpe >= 70 ) ) THEN
--> your mods
ENDIF


But I the WRF grid is not standard. How do I know the equivalent coordinates of 120-121.5E and 12-15N in terms of the i and j?

Sincerely,
Lyndz


kwerner said:
Hi Lyndz,
It looks like you're on the right track, and just need to make modifications to specify the exact area for the change. Take a look at this forum topic:
https://forum.mmm.ucar.edu/phpBB3/viewtopic.php?f=45&t=8413&p=14559#p14559
See my response from Sept. 26th, providing an example "if" test for something similar. The difference is that in the other topic, they are modifying the WRF code, as opposed to WPS, and therefore the naming convention for i,j will likely be a bit different. The example shows mods for distributed memory code, hence the "p" (for 'patch') in "ips/jps." Let me know if this helps.
 
You can use something like ncview to look at the geo_em* file, and under 2D variables you should find XLAT_M and XLONG_M. That is the lat/lon at the mass point of each grid (i.e., the center of each grid cell). As you scroll over it, it will tell you the corresponding i,j point.
 
If your WRF domain uses a map projection in which lines of constant latitude don't follow lines of constant j-index and lines of constant longitude don't follow lines of constant i-index, it might get tedious to find all (i,j) points that lie within a box that is bounded by constant longitudes and latitudes. In this case, having access to the full XLAT_M and XLONG_M arrays in your code to modify the topography might be helpful. Having only taken a quick look at the read_wrf_nc code, reading and passing the XLAT_M and XLONG_M arrays to the USER_CODE function might be a bit tedious. If you happen to be comfortable with Python, something like the following might be helpful:
Code:
from netCDF4 import Dataset
  
f = Dataset('geo_em.d01.nc', mode='r+')

xlat = f.variables['XLAT_M'][:]
xlon = f.variables['XLONG_M'][:]
hgt = f.variables['HGT_M'][:]
e_we = len(f.dimensions['west_east'])
e_sn = len(f.dimensions['south_north'])

for j in range(e_sn):
        for i in range(e_we):
                if xlat[0,j,i] >= 35.0 and xlat[0,j,i] <= 37.0 and \
                   xlon[0,j,i] >= -83.0 and xlon[0,j,i] <= -81.0:
                        hgt[0,j,i] = 0.5 * hgt[0,j,i]

f.variables['HGT_M'][:] = hgt
f.close()
The above code relies on the Python netCDF4 module. It's not the most efficient or compact code, but it's hopefully understandable enough to demonstrate the basic approach.
 
Hi @mgduda,

I was wondering too, using this method. How can I estimate the exact boundary 120-121.5E and 12-15N.
I tried to use ncview but using this method, I can only get the approximate location of the box not the exact one.

Is there a way to determine this if i used python?

Sincerely,
Lyndz
 
I may not be fully understanding the question. Since the example Python code is looping over all grid cells in the geo_em.d01.nc file and accessing the latitude and longitude coordinates for each of those grid cells, all that you need to do in each iteration of the loop is to test whether those geographical coordinates lie within you region of interest. For a region that is bounded by constant values of latitude and longitude, this is a simple comparison of the grid cell coordinates with the bounding latitude and longitude as in the example script. In general, if the region for which you would like to modify the topography is described as a polygon, you'll need to implement a method for testing whether a given latitude and longitude is within this polygon.
 
@mgduda,

Thank you for the help on this.

Do you know how I can edit this python script if I want to use a shapefile?

I uploaded my files here:
https://www.dropbox.com/sh/9cby7kegv8c35d6/AADLnbtMz5Y3Tqg-Ck2KYMywa?dl=0

Code:
from netCDF4 import Dataset
  
f = Dataset('geo_em.d01.nc', mode='r+')

xlat = f.variables['XLAT_M'][:]
xlon = f.variables['XLONG_M'][:]
hgt = f.variables['HGT_M'][:]
e_we = len(f.dimensions['west_east'])
e_sn = len(f.dimensions['south_north'])

for j in range(e_sn):
        for i in range(e_we):
                if xlat[0,j,i] >= 35.0 and xlat[0,j,i] <= 37.0 and \
                   xlon[0,j,i] >= -83.0 and xlon[0,j,i] <= -81.0:
                        hgt[0,j,i] = 0.5 * hgt[0,j,i]

f.variables['HGT_M'][:] = hgt
f.close()

Apologies. I am very new to python. This previous python script works well.
 
Top