Re: coordinate problem

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Thu Jan 12 2012 - 10:35:59 MST

Hi JinQ,

You are trying to subscript the variable with:

> u_planeb=wgt_areaave(u_plane1({a:b},{c:d}),1.,1.,0)-wgt_areaave(u_plane0({a:b},{c:d}),1.,1.,0)

You are dealing with WRF data, which doesn't have coordinate arrays. This is required in order to subscript with coordinate values.

Instead of using this, you can look at the wrf_user_ll_to_ij function. This returns an i,j index
based on a lat/lon value. You can then use the i,j to subscript your array.

--Mary

On Jan 11, 2012, at 8:18 PM, JinQ Liu wrote:

> Hi,
> I have add the "degrees_north" attributes on the variable, but still have this error :
> fatal:Dimension (west_east) of (u_plane1) does not have an associated coordinate variable
>
> Any suggestions will be appreciate.Thanks.
> JinQ
>
> --
> JinQ, M.S.
> Ocean University of China
> --
>
> ftp> put wrfoutd01_origin
>
> scripts:
> 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/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" ; WRF_Times2Udunits_c
> ;**************************************************
> begin
> xp = new((/3,6/),float) ; allocate memory
> yp = new((/3,6/),float)
> yp(0,:) = (/ 23.2, 24.3, 25.6, 26.8, 27.8, 29.0 /)
> xp(0,:) = (/ 119.0,118.7,118.7,118.7,119.2,119.9 /)
> yp(1,:) = (/ 23.382, 24.1339, 25.2679, 26.1043, 26.8851, 27.8487/)
> xp(1,:) = (/ 119., 118.945, 118.838, 118.786, 118.884, 119.472/)
> yp(2,:) = (/ 23.382, 24.1336, 25.402, 26.4258, 27.421, 28.3855/)
> xp(2,:) = (/ 119. , 118.975, 118.929, 119.03, 119.192, 119.571/)
>
> res=True
> f = addfile("wrfoutd01_origin.nc","r")
> loc=wrf_user_ll_to_ij(f,xp(1,:),yp(1,:),res)
> print(loc)& nbsp;
> lon1=loc(0,:)+11.11
> lon2=loc(0,:)-11.11
> lat1=loc(1,:)+11.11
> lat2=loc(1,:)-11.11
>
> afterloc1=wrf_user_ij_to_ll(f,lon1,lat1,res)
> afterloc2=wrf_user_ij_to_ll(f,lon1,lat2,res)
> afterloc3=wrf_user_ij_to_ll(f,lon2,lat1,res)
> afterloc4=wrf_user_ij_to_ll(f,lon2,lat2,res)
>
> wks = gsn_open_wks("x11" ,"radius")
> write_matrix(afterloc1,"6f10.2",False)
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;
> a=afterloc2(1,0)
> b=afterloc1(1,0)
> c=afterloc2(0,0)
> d=afterloc1(0,0)
> print(a)
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> times = wrf_user_list_times(f) ; get times in the file
> pressure_levels = (/ 850., 200./) ; pressure levels to plot
> nlevels = dimsizes(pressure_levels) ; number of pressure level s
> do it=3,4
> print("Working on time: " + times(it) )
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> u = wrf_user_getvar(f,"ua",it) ; u averaged to mass points
> v = wrf_user_getvar(f,"va",it) ; v averaged to mass points
> p = wrf_user_getvar(f, "pressure",it) ; pressure is our vertical coordinate
> z = wrf_user_getvar(f, "z",it) ; grid point height
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> pressure0 = pressure_levels(0)
> pressure1 = pressure_levels(1)
> u_plane0 = wrf_user_intrp3d( u,p,"h",pressure0,0.,False)
> v_plane0 ; = wrf_user_intrp3d( v,p,"h",pressure0,0.,False)
> u_plane1 = wrf_user_intrp3d( u,p,"h",pressure1,0.,False)
> v_plane1 = wrf_user_intrp3d( v,p,"h",pressure1,0.,False)
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> lat2d=f->XLAT(0,:,:)
> lon2d=f->XLONG(0,:,:)
> lat2d@units="degrees_north"
> lon2d@units="degrees_east"
> u_plane1@lat2d=lat2d
> u_plane1@lon2d=lon2d
> v_plane1@lat2d=lat2d
> v_plane1@lon2d=lon2d
> u_plane0@lat2d=lat2d
> u_plane0@lon2d=lon2d
> v_plane0@lat2d=lat2d
> v_plane0@lon2d=lon2d
>
> printVarSummary(u_plane1)
> u_planeb=wgt_areaave(u_plane1({a:b},{c:d}),1.,1.,0)-wgt_areaave(u_plane0({a:b},{c:d}),1.,1.,0)
> v_planeb=wgt_areaave(v_plane1({a:b},{c:d}),1.,1.,0)-wgt_areaave(v_plane0({a:b},{c:d}),1.,1.,0)
> u_plane = u_plane*1.94386 ; kts
> v_plane = v_plane*1.94386 ; kts
> u_plane@units = "kts"
> v_plane@units = "kts"
> print(u_plane)
> wind=((u_plane)^2+(v_plane)^2)^0.5
> print(wind)
>
> end do
> end
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jan 12 10:36:12 2012

This archive was generated by hypermail 2.1.8 : Wed Jan 18 2012 - 09:21:55 MST