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

NCL Scripts

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.

sainani7

New member
Dear all am running plotfmt.ncl script using CCSM4 data input but its shows "fatal:(lon) is not a named dimension in variable (slab).
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 199 in file plotfmt.ncl"


my file is.
;_______________________________________________________________________________________
;To run the script type:
; ncl plotfmt.ncl {input file}
;
; e.g.
; ncl plotfmt.ncl 'filename="FILE:2005-06-01_00"'
;
;
; This script can only be used in NCL V6.2 or later!!!!!
;_______________________________________________________________________________________

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

; Make sure we have a datafile to work with
if (.not. isvar("filename") ) then
print(" ")
print(" ")
print(" ")
print(" ")
print(" ")
print(" ")
print(" ")
exit
end if

head_real = new(14,float)
field = new(1,string)
hdate = new(1,string)
units = new(1,string)
map_source = new(1,string)
desc = new(1,string)

; We generate plots, but what kind do we prefer?
type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"

outf = "plotfmt_bin"
wks = gsn_open_wks(type,outf)

res = True
res@cnFillOn = True

;Open binary file

istatus = wrf_wps_open_int(filename)

if(istatus.eq.0) then

;Read header of file
wrf_wps_rdhead_int(istatus,head_real,field,hdate, \
units,map_source,desc)

do while (istatus.eq.0)

version = toint(head_real(0))
xfcst = 0
xlvl = 850
nx = toint(5)
ny = toint(5)
iproj = 1
; print("==================================================")
print("VAR = " + field + "__" + xlvl )

lat0 = 5
lon0 = 100

print("hdate ='2015-01-01_00:00:0000000'" + hdate + "'")
print("units ='" + units)
print("desc = '" + desc + "'")
print("field = '" + field + "'")
print("map_source = '" + map_source)
print("version = " + version)
print("xfcst = " + xfcst)
print("xlvl = '" + xlvl)
print("nx/ny = " + nx + "/" + ny )
print("iproj = '" + iproj)
print("lat0/lon0 = " + lat0 + "/" + lon0)
print("head_real(8) = " + head_real(8))
print("head_real(9) = " + head_real(9))
print("head_real(10) = " + head_real(10))
print("head_real(11) = " + head_real(11))
print("head_real(12) = " + head_real(12))
print("head_real(13) = " + head_real(13))


if (iproj.eq.0) then ;Cylindrical Equidistant
lat = lat0 + ispan(0,ny-1,1)*head_real(8)
lon = lon0 + ispan(0,nx-1,1)*head_real(9)
;---Turn these into 1D coordinate arrays
lat!0 = "lat"
lon!0 = "lon"
lon@units = "degrees_east"
lat@units = "degrees_north"
lat&lat = lat
lon&lon = lon
end if

if (iproj.eq.1) then ; Mercator
dx = 3
dy = 3
truelat1 = 5
xlonc = 60
res1 = True
res1@MAP_PROJ = 1
res1@TRUELAT1 = 5
res1@TRUELAT2 = 25
res1@STAND_LON = xlonc
res1@DX = dx*1000.
res1@DY = dy*1000.
res1@REF_LAT = lat0
res1@REF_LON = lon0
res1@POLE_LAT = 90.0
res1@POLE_LON = 0.0
res1@LATINC = 3
res1@LONINC = 3
res1@KNOWNI = 1.0
res1@KNOWNJ = 1.0
loc = wrf_ij_to_ll (10.0,20.0,res1)

res@gsnAddCyclic = False
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat0
res@mpLeftCornerLonF = lon0
res@mpRightCornerLatF = loc(1)
res@mpRightCornerLonF = loc(0)
res@tfDoNDCOverlay = True
res@mpProjection = "mercator"
end if

if (iproj.eq.3) then ; Lambert Conformal
dx = head_real(8)
dy = head_real(9)
truelat1 = 5
truelat2 = 25
xlonc = 60
res1 = True
res1@MAP_PROJ = 1
res1@TRUELAT1 = 5
res1@TRUELAT2 = 25
res1@STAND_LON = xlonc
res1@DX = dx*1000.
res1@DY = dy*1000.
res1@REF_LAT = lat0
res1@REF_LON = lon0
res1@POLE_LAT = 90.0
res1@POLE_LON = 0.0
res1@LATINC = 0.0
res1@LONINC = 0.0
res1@KNOWNI = 1.0
res1@KNOWNJ = 1.0
loc = wrf_ij_to_ll (nx,ny,res1)

res@gsnAddCyclic = False
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat0
res@mpLeftCornerLonF = lon0
res@mpRightCornerLatF = loc(1)
res@mpRightCornerLonF = loc(0)
res@tfDoNDCOverlay = True
res@mpProjection = "LambertConformal"
res@mpLambertParallel1F = truelat1
res@mpLambertParallel2F = truelat2
res@mpLambertMeridianF = xlonc
end if

if (iproj.eq.4) then ;Gaussian
nlats = head_real(8)
deltalon = head_real(9)
deltalat = 2.*(lat0)/(2.*nlats-1)
if (lat0 .ge. 80.) then
deltalat = -1.0*deltalat
end if
lat = lat0 + ispan(0,ny-1,1)*deltalat
lon = lon0 + ispan(0,nx-1,1)*deltalon
;---Turn these into 1D coordinate arrays
lat!0 = "lat"
lon!0 = "lon"
lon@units = "degrees_east"
lat@units = "degrees_north"
lat&lat = lat
lon&lon = lon
end if

istatus = 0

; Read 2D data

slab = wrf_wps_rddata_int(istatus,5,5)

slab@_FillValue = -1e+30
slab!1 = "80"
slab!0 = "30"
slab&lon = "60"
slab&lat = "25"

slab@units = units

slab@description = xlvl +" "+ desc
; printVarSummary(slab)

map = gsn_csm_contour_map(wks,slab,res)
delete(lat)
; delete(lon)
delete(slab)

wrf_wps_rdhead_int(istatus,head_real,field,hdate, \
units,map_source,desc)

end do

end if

end





actually am trying to draw a plot with fixed latitudes and longitudes.
anyone give me the answer about this query.

Advance Thank you!
 
Hi,
Can you attach the following:

1) namelist.wps file
2) your ungribbed intermediate file that you're using to run this

and please let me know which version of WPS you are using.

Thanks,
Kelly
 
Top