Hi,
I am new to writing NCL scripts for WRF. I'm trying to zoom into a section of my inner most grid (grid 3). With no zoom, the contours, data, basemap (country outline), and Lat/Long all see to be in correct position. However, when the zoom script is used the basemap and lat/long are no longer aligned with the data and contours. I am unsure of which bit i need to modify. Please excuse the messyness of the script.
Thanks,
Andy
;******************************************************
; Script to produce wind field output for CSIP data
;******************************************************
;andrew_barkwith_at_yahoo.co.uk
;****USER**** = user defined variable
load "WRFOptions.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "gsn_code.ncl"
load "WRFPlot.ncl"
load "WRFUserARW.ncl"
;load "SkewTFunc.ncl"
;******************************************************
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
a = addfile("../WRFV2_new/test/em_real/wrfout_d03_2005-07-19_00:00:00.nc","r")
;******************************************************
;Output Type
type = "pdf"
;******************************************************
;Output File Name
wks = gsn_open_wks(type,"wrf_CSIP")
; Debug information.
debug = False
; debug = True
icount = 0
;******************************************************
;Title
res_at_MainTitle = "WRF CSIP"
;******************************************************
; Number of time steps
FirstTime = True
times = wrf_user_list_times(a) ; retreive times in the file
ntimes = dimsizes(times) ; number of times in the file
;******************************************************
; The specific plane we want to plot data on
plane = (/ 250., 250./) ; (x,y) point for vertical plane ;*****USER*****
angle = 0.0 ; angle of cross section orientation ;*****USER*****
aspect_ratio = 1 ; Aspect ratio of output ;*****USER******
;******************************************************
; This is the big loop over all of the time periods to process.
do it = 0,ntimes-1,24 ; last number is the loop number ******USER*****
time = it
res_at_TimeLabel = times(it)
res_at_AspectRatio = aspect_ratio
res_at_PlotOrientation = "(" + plane(0)+","+plane(1) + ")" + " angle " + angle
;****************************************************
;Map Background
if (FirstTime) then
times_sav = times(it)
end if
res_at_TimeLabel = times(it)
map = wrf_map(wks,a,True)
;****************************************************
; Variables
u = wrf_user_getvar(a,"ua",time) ; ua is u averaged to mass points
v = wrf_user_getvar(a,"va",time) ; va is v averaged to mass points
w = wrf_user_getvar(a,"wa",time) ; w field
z = wrf_user_getvar(a, "Z",time) ; grid point height
ter = wrf_user_getvar(a,"HGT",time) ; need terrain height sometimes
u_plane = wrf_user_intrp3d( u,z,ter,"v",plane,angle)
v_plane = wrf_user_intrp3d( v,z,ter,"v",plane,angle)
w_plane = wrf_user_intrp3d( w,z,ter,"v",plane,angle)
;*****************************************************
;1. Cross Sections for u, v, and w
;*****************************************************
; Vertical Velocity
opts_w = res
opts_w_at_FieldTitle = w_at_description
opts_w_at_cnFillOn = True
opts_w_at_gsnSpreadColorEnd = -10
opts_w_at_tiXAxisString = "Number of Grid Points"
opts_w_at_tiYAxisString = "Height (km)"
opts_w_at_tmYLMode = "Explicit"
opts_w_at_tmYLValues = ispan(0,100,10)
opts_w_at_tmYLLabels = ispan(0,20,2)
contour_w = wrf_contour(a,wks, w_plane,opts_w)
print_opts("opts_w", opts_w, debug)
;plot
wrf_overlay(wks,contour_w,False)
print_header(icount,debug)
;*****************************************************
;u velocity
opts_u = res
opts_u_at_FieldTitle = u_at_description
opts_u_at_cnFillOn = True
opts_u_at_gsnSpreadColorEnd = -10
contour_u = wrf_contour(a,wks, u_plane,opts_u)
print_opts("opts_u", opts_u, debug)
;plot
wrf_overlay(wks,contour_u,False)
print_header(icount,debug)
;*****************************************************
;v velocity
opts_v = res
opts_v_at_FieldTitle = v_at_description
opts_v_at_cnFillOn = True
opts_v_at_gsnSpreadColorEnd = -10
contour_v = wrf_contour(a,wks, v_plane,opts_v)
print_opts("opts_v", opts_v, debug)
;plot
wrf_overlay(wks,contour_v,False)
print_header(icount,debug)
;*****************************************************
;Delete some used variables
delete(u_plane)
delete(v_plane)
;*****************************************************
;2. Map View of z, u and v at heights 300, 500, 1000m
;*****************************************************
height_levels = (/ 300., 500., 1000./) ; heigth levels to plot(m) - remember to change levels further down *****USER*****
nlevels = dimsizes(height_levels) ; number of height levels
do level = 0,nlevels-1
height = height_levels(level)
u_plane = wrf_user_intrp3d( u,z,ter,"h",height,0.)
v_plane = wrf_user_intrp3d( v,z,ter,"h",height,0.)
;*****************************************************
;Terrain contours
opts_ter = res
; opts_ter_at_cnFillOn = True ; contour plot of terrain
contour_ter = wrf_contour(a,wks,ter,opts_ter)
;******************************************************
;Zoom
zooom = 1 ;zoom on=1 off=0 ;*****USER*****
if ( zooom .eq. 1)
dims = dimsizes(ter)
; As an example, we are looking for the lower right 1/4 of the domain
x_start = dims(1)/2
x_end = dims(1)-1
y_start = 0
y_end = dims(0)/2
ter_zoom = ter(y_start:y_end,x_start:x_end)
ter_zoom_at_description = ter_at_description
contour_ter = wrf_contour(a,wks,ter_zoom,opts_ter)
contour_ter_at_description = ter_at_description
u_plane_zoom = u_plane(y_start:y_end,x_start:x_end)
u_plane_zoom_at_description = u_plane_at_description
v_plane_zoom = v_plane(y_start:y_end,x_start:x_end)
v_plane_zoom_at_description = v_plane_at_description
map_zoom = wrf_map_zoom(wks,a,res,y_start,y_end,x_start,x_end)
end if
;**************************************************************
;Wind u and v as contour plots
;u
opts_windu = res
opts_windu_at_FieldTitle = "Winds"
opts_windu_at_PlotLevelID = 0.001*height + " km"
opts_windu_at_FieldTitle = u_at_description
opts_windu_at_cnFillOn = True
opts_windu_at_gsnSpreadColorEnd = -10
contour_windu = wrf_contour(a,wks,u_plane,opts_windu)
print_opts("opts_windu", opts_windu, debug)
;v
opts_windv = res
opts_windv_at_FieldTitle = "Winds"
opts_windv_at_PlotLevelID = 0.001*height + " km"
opts_windv_at_FieldTitle = v_at_description
opts_windv_at_cnFillOn = True
opts_windv_at_gsnSpreadColorEnd = -10
contour_windv = wrf_contour(a,wks,v_plane,opts_windv)
print_opts("opts_windv", opts_windv, debug)
;*******************************************************
;300m
if ( height .eq. 300 ) then
if ( zooom .eq. 1) then ; zoom on
;u
contour_windu = wrf_contour(a,wks,u_plane_zoom,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windu/),True)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane_zoom,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windv/),True)
print_header(icount,debug)
end if
if (zooom .eq. 0) then ; zoom off
;u
contour_windu = wrf_contour(a,wks,u_plane,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map,(/contour_ter,contour_windu/),False)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map,(/contour_ter,contour_windv/),False)
print_header(icount,debug)
end if
end if
;********************************************************
;500m
if ( height .eq. 500 ) then
if ( zooom .eq. 1) then ; zoom on
;u
contour_windu = wrf_contour(a,wks,u_plane_zoom,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windu/),True)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane_zoom,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windv/),True)
print_header(icount,debug)
end if
if (zooom .eq. 0) then ; zoom off
;u
contour_windu = wrf_contour(a,wks,u_plane,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map,contour_windu,False)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map,contour_windv,False)
print_header(icount,debug)
end if
end if
;******************************************************
;1000m
if ( height .eq. 1000 ) then
if ( zooom .eq. 1) then ;zoom on
;u
contour_windu = wrf_contour(a,wks,u_plane_zoom,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windu/),True)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane_zoom,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map_zoom,(/contour_ter,contour_windv/),True)
print_header(icount,debug)
end if
if (zooom .eq. 0) then ;zoom off
;u
contour_windu = wrf_contour(a,wks,u_plane,opts_windu)
print_opts("opts_windu", opts_windu, debug)
wrf_map_overlay(wks,map,contour_windu,False)
print_header(icount,debug)
;v
contour_windv = wrf_contour(a,wks,v_plane,opts_windv)
print_opts("opts_windv", opts_windv, debug)
wrf_map_overlay(wks,map,contour_windv,False)
print_header(icount,debug)
end if
end if
end do
end do
end
__________________________________________________________________________________________________________
Andrew Barkwith
(MEarthSci)
0780 264 9917
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 24 2008 - 10:51:51 MST
This archive was generated by hypermail 2.2.0 : Tue Nov 25 2008 - 10:18:44 MST