A dip into NCL

NCL is an interpreted language designed for scientific data analysis and visualization. I use it mostly for visualization, but not seriously, because there is another developed platform for visualization. For NCL, I just use it to check the model results.

So I just got a dip of NCL for the last three days, in an attempt to plot a wind vertical profile with eta level. The wind vector goes along the terrain and the plot is blank where there is a hill, so that the impact of topography could be seen. Shame on myself, I still haven't gone through it. And in terms of time, I have to give up. Here I'd lie to summarize up what I've learned with NCL.

Coordinate variables are the key information to plot right. Variables from model results (in my case, WRF, CMAQ, and SMOKE) usually bear no coordinate, but columns and rows. It is fine if you just want to see the patterns. However, if you want a base map overlaid, the right coordinates have to be assigned. Coordinates are usually stored in another output variable, e.g. XLAT for latitude and XLONG for longitude in WRF. What you have to do is to assign lat/lon to the right dimension of, say, wind. Here is the code:

f         = addfile(filename, "r")
U        = f->U                             ; wind in east-west direction, time*lev*lat*lon
W        = f->W                            ; wind in bottom-up direction
lat       = f->XLAT(0,:,0)             ; time*lat*lon
lon      = f->XLON(0,0,:)
znu     = f->ZNU(0,:)                 ; time*lev

lat@units  = "degrees_north" 
lon@units = "degrees_east"
lat!0          = "lat"
lon!0         = "lon"
lat&lat      =  lat
lon&lon     =  lon

lev                      = znu*1000            ; [-105.1526..-82.84741]
lev@long_name  = "eta*1000"
lev@units           = "hPa"
lev!0                   = "lev"
lev&lev               =  lev

U!0      = "lev"
U!1      = "lat"
U!2      = "lon"
U&lev    =  lev
U&lat    =  lat
U&lon    =  lon

; And the same for w, left out here

But wait, there is a mistake here. NCL would prompt dimension inconsistent between U and W. Let's take a close look at it. By the funtion printVarSummary we find the dimension of both U and W is time*level*latitude*longitude. So, why inconsistent?

It comes out that there are two grids in WRF, the staggered grid and the mass grid. Value in mass grid is refers to value in the center of the grid while value in staggered grid means value in the boundary of the grid. So the dimension size of staggered grids is always greater than mass grids by 1. To solve the inconsistency, we average the value to mass grid points.

dimU = dimsizes(U)
nlonU = dimU(3)
u = 0.5 * (U(:,:,:,0:nlonU-2) + U(:,:,:,1:nlonU-1))

In this case, u is on the common grid. And the same is for w.

Finally gsn_csm_pres_hgt_vector (example is here and here) is used to plot the vector. However, the result (shown in figure below, at 24.5N along 113-114E) is quite different from the example, mainly for the three points:

 1. I am not sure whether the eta level is terrain-following. Though there are "ups and downs", no specific "hills" are to be found. Of course, it may be a problem of scales.
2. The temperature contour goes along eta level, which is weird.
3. Wind always flows to the east. There seems no "up" wind.

Some other resources:
To plot vectors:
gsn_vector(), gsn_csm_vector(). And the latter is more advanced.
gsn_csm_vector_map() overlays the vector plot on a base map.

To draw vector and scalar simultaneously.
gsn_csm_vector_scalar_map()