From 0983f65853b277b272ba9ce10698cbdb1322c575 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Sat, 22 Mar 2014 00:29:06 +0000 Subject: [PATCH] Add an NCL script for plotting intermediate files. This script requires NCL 6.2 or later. Thanks to Abby Jaye. A util/plotfmt.ncl git-svn-id: https://svn-wrf-wps.cgd.ucar.edu/trunk@811 86b71a92-4018-0410-97f8-d555beccfc3a --- util/plotfmt.ncl | 215 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 util/plotfmt.ncl diff --git a/util/plotfmt.ncl b/util/plotfmt.ncl new file mode 100644 index 0000000..8502cdc --- /dev/null +++ b/util/plotfmt.ncl @@ -0,0 +1,215 @@ +;_______________________________________________________________________________________ +;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(" ### MUST SUPPLY a filename ### ") + print(" ") + print(" Something like: ") + print(" ncl plotfmt.ncl filename=FILE:2005-06-01_00") + print(" REMEMBER TO ADD QUOTES" ) + print(" Refer to the information at the top of this file for more info and syntax" ) + 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 = head_real(1) + xlvl = head_real(2) + nx = toint(head_real(3)) + ny = toint(head_real(4)) + iproj = toint(head_real(5)) + print("==================================================") + print("VAR = " + field + "__" + xlvl ) + + lat0 = head_real(6) + lon0 = head_real(7) + + print("hdate = '" + 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 = head_real(8) + dy = head_real(9) + truelat1 = head_real(10) + res1 = True + res1@MAP_PROJ = 3 + res1@TRUELAT1 = truelat1 + 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 = "mercator" + end if + + if (iproj.eq.3) then ; Lambert Conformal + dx = head_real(8) + dy = head_real(9) + xlonc = head_real(10) + truelat1 = head_real(11) + truelat2 = head_real(12) + res1 = True + res1@MAP_PROJ = 1 + res1@TRUELAT1 = truelat1 + res1@TRUELAT2 = truelat2 + 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,nx,ny) + + slab@_FillValue = -1e+30 + slab!1 = "lon" + slab!0 = "lat" + slab&lon = lon + slab&lat = lat + 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 -- 2.11.4.GIT