Hi Andrew,
I have an idea of what might be the problem, but I need your script and
data. Can you email me separately and let me know if you can provide
these files? I will provide an ftp location if you need one.
--Mary
On Mon, 24 Nov 2008, Andrew Barkwith wrote:
> 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 - 14:44:49 MST
This archive was generated by hypermail 2.2.0 : Tue Nov 25 2008 - 10:18:44 MST