util/plotfmt.ncl

   1 ;_______________________________________________________________________________________
   2 ;To run the script type:
   3 ;       ncl plotfmt.ncl {input file}
   4 ;
   5 ;       e.g.
   6 ;               ncl plotfmt.ncl 'filename="FILE:2005-06-01_00"'
   7 ;
   8 ;
   9 ;  This script can only be used in NCL V6.2 or later!!!!!
  10 ;_______________________________________________________________________________________
  11
  12 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  13 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  14 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
  15
  16 begin
  17
  18 ; Make sure we have a datafile to work with
  19   if (.not. isvar("filename") ) then
  20     print(" ")
  21     print(" ### MUST SUPPLY a filename ### ")
  22     print(" ")
  23     print("     Something like: ")
  24     print("     ncl plotfmt.ncl filename=FILE:2005-06-01_00")
  25     print("          REMEMBER TO ADD QUOTES" )
  26     print("          Refer to the information at the top of this file for more info and syntax" )
  27     exit
  28   end if
  29
  30   head_real  = new(14,float)
  31   field      = new(1,string)
  32   hdate      = new(1,string)
  33   units      = new(1,string)
  34   map_source = new(1,string)
  35   desc       = new(1,string)
  36
  37 ; We generate plots, but what kind do we prefer?
  38  type = "x11"
  39 ; type = "pdf"
  40 ; type = "ps"
  41 ; type = "ncgm"
  42
  43   outf = "plotfmt_bin"
  44   wks = gsn_open_wks(type,outf)
  45
  46   res = True
  47   res@cnFillOn = True
  48
  49   ;Open binary file
  50
  51   istatus =  wrf_wps_open_int(filename)
  52
  53   if(istatus.eq.0) then
  54
  55     ;Read header of file
  56    wrf_wps_rdhead_int(istatus,head_real,field,hdate, \
  57                       units,map_source,desc)
  58
  59     do while (istatus.eq.0)
  60
  61       version = toint(head_real(0))
  62       xfcst   = head_real(1)
  63       xlvl    = head_real(2)
  64       nx      = toint(head_real(3))
  65       ny      = toint(head_real(4))
  66       iproj   = toint(head_real(5))
  67       print("==================================================")
  68       print("VAR = " + field + "__" + xlvl )
  69
  70       lat0      = head_real(6)
  71       lon0      = head_real(7)
  72
  73       print("hdate         = '" + hdate + "'")
  74       print("units         = '" + units + "'")
  75       print("desc          = '" + desc + "'")
  76       print("field         = '" + field + "'")
  77       print("map_source    = '" + map_source + "'")
  78       print("version       = " + version)
  79       print("xfcst         = " + xfcst)
  80       print("xlvl          = " + xlvl)
  81       print("nx/ny         = " + nx + "/" + ny )
  82       print("iproj         = " + iproj)
  83       print("lat0/lon0     = " + lat0 + "/" + lon0)
  84       print("head_real(8)  = " + head_real(8))
  85       print("head_real(9)  = " + head_real(9))
  86       print("head_real(10) = " + head_real(10))
  87       print("head_real(11) = " + head_real(11))
  88       print("head_real(12) = " + head_real(12))
  89       print("head_real(13) = " + head_real(13))
  90
  91
  92       if (iproj.eq.0) then         ;Cylindrical Equidistant
  93         lat = lat0 + ispan(0,ny-1,1)*head_real(8)
  94         lon = lon0 + ispan(0,nx-1,1)*head_real(9)
  95 ;---Turn these into 1D coordinate arrays
  96         lat!0     = "lat"
  97         lon!0     = "lon"
  98         lon@units = "degrees_east"
  99         lat@units = "degrees_north"
 100         lat&lat   = lat
 101         lon&lon   = lon
 102       end if
 103
 104       if (iproj.eq.1) then          ; Mercator
 105         dx = head_real(8)
 106         dy = head_real(9)
 107         truelat1 = head_real(10)
 108         res1 = True
 109         res1@MAP_PROJ  = 3
 110         res1@TRUELAT1  = truelat1
 111         res1@DX        = dx*1000.
 112         res1@DY        = dy*1000.
 113         res1@REF_LAT   = lat0
 114         res1@REF_LON   = lon0
 115         res1@POLE_LAT  = 90.0
 116         res1@POLE_LON  =  0.0
 117         res1@LATINC    = 0.0
 118         res1@LONINC    = 0.0
 119         res1@KNOWNI    = 1.0
 120         res1@KNOWNJ    = 1.0
 121         loc = wrf_ij_to_ll (nx,ny,res1)
 122
 123         res@gsnAddCyclic = False
 124         res@mpLimitMode = "Corners"
 125         res@mpLeftCornerLatF = lat0
 126         res@mpLeftCornerLonF = lon0
 127         res@mpRightCornerLatF = loc(1)
 128         res@mpRightCornerLonF = loc(0)
 129         res@tfDoNDCOverlay = True
 130         res@mpProjection = "mercator"
 131       end if
 132
 133       if (iproj.eq.3) then          ; Lambert Conformal
 134         dx = head_real(8)
 135         dy = head_real(9)
 136         xlonc = head_real(10)
 137         truelat1 = head_real(11)
 138         truelat2 = head_real(12)
 139         res1 = True
 140         res1@MAP_PROJ  = 1
 141         res1@TRUELAT1  = truelat1
 142         res1@TRUELAT2  = truelat2
 143         res1@STAND_LON = xlonc
 144         res1@DX        = dx*1000.
 145         res1@DY        = dy*1000.
 146         res1@REF_LAT   = lat0
 147         res1@REF_LON   = lon0
 148         res1@POLE_LAT  = 90.0
 149         res1@POLE_LON  =  0.0
 150         res1@LATINC    = 0.0
 151         res1@LONINC    = 0.0
 152         res1@KNOWNI    = 1.0
 153         res1@KNOWNJ    = 1.0
 154         loc = wrf_ij_to_ll (nx,ny,res1)
 155
 156         res@gsnAddCyclic = False
 157         res@mpLimitMode = "Corners"
 158         res@mpLeftCornerLatF = lat0
 159         res@mpLeftCornerLonF = lon0
 160         res@mpRightCornerLatF = loc(1)
 161         res@mpRightCornerLonF = loc(0)
 162         res@tfDoNDCOverlay = True
 163         res@mpProjection = "LambertConformal"
 164         res@mpLambertParallel1F = truelat1
 165         res@mpLambertParallel2F = truelat2
 166         res@mpLambertMeridianF = xlonc
 167         res@pmTickMarkDisplayMode = "Always"
 168         res@mpGridAndLimbOn = True
 169       end if
 170
 171       if (iproj.eq.4) then        ;Gaussian
 172         nlats    = head_real(8)
 173         deltalon = head_real(9)
 174         deltalat = 2.*(lat0)/(2.*nlats-1)
 175         if (lat0 .ge. 80.) then
 176           deltalat = -1.0*deltalat
 177         end if
 178         lat = lat0 + ispan(0,ny-1,1)*deltalat
 179         lon = lon0 + ispan(0,nx-1,1)*deltalon
 180 ;---Turn these into 1D coordinate arrays
 181         lat!0     = "lat"
 182         lon!0     = "lon"
 183         lon@units = "degrees_east"
 184         lat@units = "degrees_north"
 185         lat&lat   = lat
 186         lon&lon   = lon
 187       end if
 188
 189       istatus = 0
 190
 191       ; Read 2D data
 192
 193       slab = wrf_wps_rddata_int(istatus,nx,ny)
 194
 195       slab@_FillValue = -1e+30
 196       if (iproj.ne.3) then ; NOT Lambert Conformal
 197          slab!1 = "lon"
 198          slab!0 = "lat"
 199          slab&lon = lon
 200          slab&lat = lat
 201       end if
 202       slab@units = units
 203
 204       slab@description = xlvl +"  "+ desc
 205 ;     printVarSummary(slab)
 206
 207       map = gsn_csm_contour_map(wks,slab,res)
 208       if (iproj.ne.3) then ; NOT Lambert Conformal
 209          delete(lat)
 210          delete(lon)
 211       end if
 212       delete(slab)
 213
 214       wrf_wps_rdhead_int(istatus,head_real,field,hdate, \
 215                          units,map_source,desc)
 216
 217     end do
 218
 219   end if
 220
 221 end