util/plotgrids_new.ncl

   1
   2 ;   Script display location of model domains
   3 ;   Only works for ARW domains
   4 ;   Only works for NCL versions later than 6.1
   5 ;   Reads namelist file directly
   6
   7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   8 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
   9
  10 begin
  11 ;
  12
  13 ; Check the version of NCL
  14   version = systemfunc("ncl -V")
  15   if(version.lt.6.2) then
  16     print("You need NCL V6.2 or later to run this script. Try running plotgrids.ncl. Stopping now...")
  17     return
  18   end if
  19
  20 ; We generate plots, but what kind do we prefer?
  21   type = "x11"
  22 ; type = "pdf"
  23 ; type = "ps"
  24 ; type = "ncgm"
  25   wks = gsn_open_wks(type,"wps_show_dom")
  26
  27 ; read the following namelist file
  28   filename = "namelist.wps"
  29
  30 ; Set the colors to be used
  31   colors = (/"white","black","White","ForestGreen","DeepSkyBlue","Red","Blue"/)
  32   gsn_define_colormap(wks, colors)
  33
  34
  35 ; Set some map information ; line and text information
  36   mpres = True
  37   mpres@mpFillOn = True
  38   mpres@mpFillColors  = (/"background","DeepSkyBlue","ForestGreen","DeepSkyBlue", "transparent"/)
  39   mpres@mpDataBaseVersion           = "Ncarg4_1"
  40   mpres@mpGeophysicalLineColor      = "Black"
  41   mpres@mpGridLineColor             = "Black"
  42   mpres@mpLimbLineColor             = "Black"
  43   mpres@mpNationalLineColor         = "Black"
  44   mpres@mpPerimLineColor            = "Black"
  45   mpres@mpUSStateLineColor          = "Black"
  46 ;  mpres@mpOutlineBoundarySets       = "AllBoundaries"
  47   ;mpres@mpGridSpacingF              = 45
  48   mpres@tiMainString                = " WPS Domain Configuration  "
  49
  50   lnres = True
  51   lnres@gsLineThicknessF = 2.5
  52   lnres@domLineColors    = (/ "white", "Red" , "Red" , "Blue" /)
  53
  54   txres = True
  55   txres@txFont = "helvetica-bold"
  56   ;txres@txJust = "BottomLeft"
  57   txres@txJust = "TopLeft"
  58   txres@txPerimOn = False
  59   txres@txFontHeightF = 0.015
  60
  61 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  62 ; Do not change anything between the ";;;;;" lines
  63
  64   maxdom = 21
  65   nvar = 19
  66   parent_idn = new (maxdom,integer)
  67   parent_grid_ration = new (maxdom,integer)
  68   i_parent_startn = new (maxdom,integer)
  69   j_parent_startn = new (maxdom,integer)
  70   e_wen = new (maxdom,integer)
  71   e_snn = new (maxdom,integer)
  72   plotvar = new((/maxdom,nvar/),float)
  73   plotvar@_FillValue = -999.0
  74
  75   plotvar = wrf_wps_read_nml(filename)
  76
  77   mpres@max_dom = floattointeger(plotvar(0,0))
  78   mpres@dx = plotvar(0,1)
  79   mpres@dy = plotvar(0,2)
  80   if (.not.ismissing(plotvar(0,3))) then
  81     mpres@ref_lat = plotvar(0,3)
  82   else
  83     mpres@ref_lat = 0.0
  84   end if
  85   if (.not.ismissing(plotvar(0,4))) then
  86     mpres@ref_lon = plotvar(0,4)
  87   else
  88     mpres@ref_lon = 0.0
  89   end if
  90   if (.not.ismissing(plotvar(0,5))) then
  91     mpres@ref_x = plotvar(0,5)
  92   end if
  93   if (.not.ismissing(plotvar(0,6))) then
  94     mpres@ref_y = plotvar(0,6)
  95   end if
  96   mpres@truelat1 = plotvar(0,7)
  97   mpres@truelat2 = plotvar(0,8)
  98   mpres@stand_lon = plotvar(0,9)
  99   mproj_int = plotvar(0,10)
 100   mpres@pole_lat = plotvar(0,11)
 101   mpres@pole_lon = plotvar(0,12)
 102
 103   do i = 0,maxdom-1
 104     parent_idn(i) = floattointeger(plotvar(i,13))
 105     parent_grid_ration(i) = floattointeger(plotvar(i,14))
 106     i_parent_startn(i) = floattointeger(plotvar(i,15))
 107     j_parent_startn(i) = floattointeger(plotvar(i,16))
 108     e_wen(i) = floattointeger(plotvar(i,17))
 109     e_snn(i) = floattointeger(plotvar(i,18))
 110   end do
 111
 112
 113   if(mpres@max_dom .gt. 1) then
 114     do i = 1,mpres@max_dom-1
 115
 116       ;Making sure edge is nested grid is at least 5 grid points from mother domain.
 117       if(i_parent_startn(i) .lt. 5) then
 118         print("Warning: Western edge of grid must be at least 5 grid points from mother domain!")
 119       end if
 120       if(j_parent_startn(i) .lt. 5) then
 121         print("Warning: Southern edge of grid must be at least 5 grid points from mother domain!")
 122       end if
 123       pointwe = (e_wen(i)-1.)/parent_grid_ration(i)
 124       pointsn = (e_snn(i)-1.)/parent_grid_ration(i)
 125       gridwe = e_wen(parent_idn(i)-1)-(pointwe+i_parent_startn(i))
 126       gridsn = e_snn(parent_idn(i)-1)-(pointsn+j_parent_startn(i))
 127       if(gridwe .lt. 5) then
 128         print("Warning: Eastern edge of grid must be at least 5 grid points from mother domain!")
 129       end if
 130       if(gridsn .lt. 5) then
 131         print("Warning: Northern edge of grid must be at least 5 grid points from mother domain!")
 132       end if
 133
 134       ;Making sure nested grid is fully contained in mother domain.
 135       gridsizewe = (((e_wen(parent_idn(i)-1)-4)-i_parent_startn(i))*parent_grid_ration(i))-(parent_grid_ration(i)-1)
 136       gridsizesn = (((e_snn(parent_idn(i)-1)-4)-j_parent_startn(i))*parent_grid_ration(i))-(parent_grid_ration(i)-1)
 137       if(gridwe .lt. 5) then
 138         print("Warning: Inner nest (domain = " + (i+1) + ") is not fully contained in mother nest (domain = " + parent_idn(i) + ")!")
 139         print("For the current setup of mother domain = " + parent_idn(i) + ", you can only have a nest of size " + gridsizewe + "X" + gridsizesn + ". Stopping Program!")
 140         exit
 141       end if
 142       if(gridsn .lt. 5) then
 143         print("Warning: Inner nest (domain = " + (i+1) + ") is not fully contained in mother nest (domain = " + parent_idn(i) + ")!")
 144         print("For the current setup of mother domain = " + parent_idn(i) + ", you can only have a nest of size " + gridsizewe + "X" + gridsizesn + ". Stopping Program!")
 145         exit
 146       end if
 147
 148       ;Making sure the nest ends at a mother grid domain point.
 149       pointwetrunc = decimalPlaces(pointwe,0,False)
 150       pointsntrunc = decimalPlaces(pointsn,0,False)
 151       if((pointwe-pointwetrunc) .ne. 0.) then
 152         nest_we_up = (ceil(pointwe)*parent_grid_ration(i))+1
 153         nest_we_dn = (floor(pointwe)*parent_grid_ration(i))+1
 154         print("Nest does not end on mother grid domain point. Try " + nest_we_dn + " or " + nest_we_up + ".")
 155       end if
 156       if((pointsn-pointsntrunc) .ne. 0.) then
 157         nest_sn_up = (ceil(pointsn)*parent_grid_ration(i))+1
 158         nest_sn_dn = (floor(pointsn)*parent_grid_ration(i))+1
 159         print("Nest does not end on mother grid domain point. Try " + nest_sn_dn + " or " + nest_sn_up + ".")
 160       end if
 161
 162     end do
 163   end if
 164
 165   mpres@parent_id = parent_idn(0:mpres@max_dom-1)
 166   mpres@parent_grid_ratio = parent_grid_ration(0:mpres@max_dom-1)
 167   mpres@i_parent_start = i_parent_startn(0:mpres@max_dom-1)
 168   mpres@j_parent_start = j_parent_startn(0:mpres@max_dom-1)
 169   mpres@e_we = e_wen(0:mpres@max_dom-1)
 170   mpres@e_sn = e_snn(0:mpres@max_dom-1)
 171
 172   if(mproj_int .eq. 1) then
 173     mpres@map_proj = "lambert"
 174     mpres@pole_lat = 0.0
 175     mpres@pole_lon = 0.0
 176   else if(mproj_int .eq. 2) then
 177     mpres@map_proj = "mercator"
 178     mpres@pole_lat = 0.0
 179     mpres@pole_lon = 0.0
 180   else if(mproj_int .eq. 3) then
 181     mpres@map_proj = "polar"
 182     mpres@pole_lat = 0.0
 183     mpres@pole_lon = 0.0
 184   else if(mproj_int .eq. 4) then
 185     mpres@map_proj = "lat-lon"
 186   end if
 187   end if
 188   end if
 189   end if
 190
 191 ; Deal with global wrf domains that don't have dx or dy
 192
 193   if (mpres@dx.lt.1e-10 .and. mpres@dx.lt.1e-10) then
 194     mpres@dx = 360./(mpres@e_we(0) - 1)
 195     mpres@dy = 180./(mpres@e_sn(0) - 1)
 196     mpres@ref_lat = 0.0
 197     mpres@ref_lon = 180.0
 198   end if
 199
 200   mp = wrf_wps_dom (wks,mpres,lnres,txres)
 201
 202 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 203
 204 ; Now you can add some information to the plot.
 205 ; Below is an example of adding a white dot over the DC location.
 206   ;pmres = True
 207   ;pmres@gsMarkerColor = "White"
 208   ;pmres@gsMarkerIndex = 16
 209   ;pmres@gsMarkerSizeF = 0.01
 210   ;gsn_polymarker(wks,mp,-77.26,38.56,pmres)
 211
 212
 213   frame(wks)           ; lets frame the plot - do not delete
 214
 215 end