Missing values using wrf_interp_3d to interpolate u & v.

The NCL graphics package.

Missing values using wrf_interp_3d to interpolate u & v.

Postby slam » Sun Jan 08, 2017 12:18 pm

Dear All,
I am trying to extract the u and v at 30m above ground levels. I tried this code in a scenario that levels were set in eta levels, the results were fine. However, I am getting missing values while running this code in scenario where model was set as pressure level.
Any thoughts? Much appreciate your help in advance.

Code: Select all
;; Load necessary libraries
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/pel/datefix.ncl"

begin

   ;; ncl doesn't seem to like reading file paths from the command line
   ;; haven't had time to work out what the fix is.
   if (.not. isvar("fpath")) then
           fpath = "wrfout_d05_2007-01-28_07:00:00"
   end if


   a = addfile(fpath+".nc","r")
   
   ;; Check for coordinates on command line
   if (.not. isvar("lat_of_interest")) then      ; is lat_of_interest on command line?
      print ("error: no latitude supplied")
      exit
   end if

   if (.not. isvar("lon_of_interest")) then      ; is lon_of_interest on command line?
        print ("error: no longitude supplied")
        exit
   end if

   ;; What times and how many time steps are in the data set?
   times  = wrf_user_list_times(a)  ; get times in the file
   ntimes = dimsizes(times)         ; number of times in the file

   wrf_lat = a->XLAT(0,:,:)
   wrf_lon = a->XLONG(0,:,:)

   diff_lat = abs(abs(lat_of_interest) - abs(wrf_lat))
   diff_lon = abs(lon_of_interest - wrf_lon)
   sum_diff = diff_lat + diff_lon

   ;; Jatin's get index from lat long script
   sum_diff_1D = ndtooned(sum_diff)
   dsizes_sum_diff = dimsizes(sum_diff)
   ;; index min provides indes of closest point
   index_min = ind_resolve(minind(sum_diff_1D),dsizes_sum_diff)

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   ; unused, but here in case they fix they're script
   ; official NCL get index from lat long script
   ;loci_ij = wrf_user_ll_to_ij(a,lon_of_interest,lat_of_interest,True)
   ;loci_ij = loci_ij -1
   ;locX = loci_ij(0)
   ;locY = loci_ij(1)
   ;print(wrf_lat(locY,locX))
   ;print(wrf_lon(locY,locX))
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

   ;; get timestep
   ts_ntimes = a->Times
   ;ts_stimes = chartostring(ts_ntimes)
   ts_dtimes = wrf_times_c(ts_ntimes, 2)
   ;timestep = convert_base10_date(ts_dtimes(1)) - convert_base10_date(ts_dtimes(0))

   ;; print headers
   print ("Time,temp_30m,rainfall,u,v,rad,mixht,sst")
   
   ;; Collect and print values
   previous_rainfall = 0.0
   do it = 0,ntimes-1,1             ; TIME LOOP

      ;temp_2m = wrf_user_getvar(a, "T2", it)
      ;temp_2m_profile = temp_2m(index_min(0,0), index_min(0,1))
               
                T_sea = wrf_user_getvar(a, "tc", it)

      mixht = wrf_user_getvar(a, "PBLH", it)
      mixht_profile = mixht(index_min(0,0), index_min(0,1))

      rainfall = wrf_user_getvar(a, "RAINNC", it)
      rainfall_profile = rainfall(index_min(0,0), index_min(0,1))
      rainfallc = wrf_user_getvar(a, "RAINC", it)
      rainfallc_profile = rainfallc(index_min(0,0), index_min(0,1))
      rainfall_total = rainfallc_profile + rainfall_profile - previous_rainfall

      rad = wrf_user_getvar(a, "SWDOWN", it)
      rad_profile = rad(index_min(0,0), index_min(0,1))

      sfcpres = wrf_user_getvar(a, "SST", it)
      sfcpres_profile = sfcpres(index_min(0,0), index_min(0,1))

                UVr = wrf_user_getvar(a, "uvmet", it)
                umet = UVr(0,:,:,:)
                vmet = UVr(1,:,:,:)
               
                ;extract at certain height
                height = wrf_user_getvar(a, "z", it)
                ter = wrf_user_getvar(a, "ter", it)
                ; Conform data to Terrain Height
                nheight = conform(height,ter,(/1,2/)) ; assuming height is a 3d array and ter is a 2d array
                height = height - nheight
               
                t_plane = wrf_user_intrp3d( T_sea,height,"h", 30,0.,False) 
                u_plane = wrf_user_intrp3d( umet,height,"h", 30,0.,False)
                v_plane  = wrf_user_intrp3d( vmet,height,"h", 30,0.,False)

                temp_30m_profile = t_plane(index_min(0,0),index_min(0,1))
                U_profiler = u_plane(index_min(0,0),index_min(0,1))
                V_profiler = v_plane(index_min(0,0),index_min(0,1))

      ;UV = wrf_user_getvar(a,"uvmet10",it)
      ;U_profile = UV(0,index_min(0,0),index_min(0,1))
      ;V_profile = UV(1,index_min(0,0),index_min(0,1)) 

      ; print out profiles
      ; print ("data: " + rainfall_profile + "," + rainfallc_profile + "," + rainfall_total + "," + previous_rainfall)
      print (times(it) + "," + temp_30m_profile + "," + rainfall_total + "," + U_profiler + "," + V_profiler + "," + rad_profile + "," + mixht_profile + "," + sfcpres_profile)
      previous_rainfall = previous_rainfall + rainfall_total

   end do

end
slam
 
Posts: 2
Joined: Sun Jan 08, 2017 3:12 am

Return to NCL

Who is online

Users browsing this forum: No registered users and 0 guests