sorting daily precipitation data.

Dennis/NCL users,

 I'm running into some sorting issues with this script. I'm working with
daily accumulated precipitation, and looking for a method to sort the data
to ultimately analyze the top 10% extreme precipitation events. I need each
precipitation data to have the associated dimensions (time,yc,xc) along with
it, to determine when these events occur.

However, output of PrcExtreme data results in all zeros (probably indicating
that it's reading the bottom 10%). Also, the time output is showing the last
10% of the entire time file, so the metadata is not being carried along with
the associated precipitation variable. I'm aware that the dim_pqsort_n
function for this script sorts both x and time, but I was wondering if
anyone had a clever method around this.

begin

;********************************************************
; Read original NetCDF file. Input necessary
; files
;*********************************************************

  fname = ("region_output.nc")
  fin = addfile (fname,"r")
  prc = fin->pr ; (time,yc, xc) 3 hrly
  lat = fin->lats_region ; lat(yc, xc)
  lon = fin->lons_region ; lon(yc, xc)
  time = fin->time
  Lambert_Conformal = fin->Lambert_Conformal

  dimx3 = dimsizes(prc)
  NTIM = dimx3(0) ; TOTAL number of time steps
  print("NTIM="+NTIM)

;************************************************************
; Create daily averages.
;************************************************************

    ntJump = 8

    x = prc(2::ntJump,:,:) ; trick: create array with meta
    printVarSummary(x) ; values will be overwritten with averages

    ntStrt = 3
    ntLast = ntStrt + ntJump - 1

    do nt=0,NTIM-9,8
       x(nt/8,:,:) = (dim_avg_n_Wrap(prc(ntStrt:ntLast,:,:), 0))
       ntStrt = ntStrt + ntJump
       ntLast = ntLast + ntJump
    end do
    printVarSummary(x)
    printMinMax(x,True)

    printMinMax(time,True)

    x@info_tag = "daily average"

;**************************************************************
; Convert units kg/(m2-s) to mm/day
;**************************************************************

    x = x*86400.
    x@units = "mm/day"
    printVarSummary(x)
    printMinMax(x,True)

;**************************************************************
; Sorting subroutine. Add if your looking for extreme
; events!
;*************************************************************

   ip = dim_pqsort_n(x, -2, 0) ; 'x' will be in descending order
   delete(ip)

   dimx = dimsizes(x)
   ntim = dimx(0)

   pc = 10 ; 10%
   prStrt = floattoint( ((100-pc)/100.)*ntim ) ; index value

;*********************************************************
; Create a netCDF with the top 10% of values.
;*********************************************************

   diro = "./"
   filo = "Precip_upper_"+pc+".nc"

   system("/bin/rm -f "+diro+filo) ; remove any pre-existing file
   ncdf = addfile(diro+filo ,"c") ; open output netCDF file

   PrcExtreme = x(prStrt::,:,:)

   ncdf->PrcExtreme = PrcExtreme ; top 'pc'% only
   ncdf->lat = lat
   ncdf->lon = lon
   ncdf->Lambert_Conformal = Lambert_Conformal

   printMinMax(PrcExtreme, True)

   globeAtt = 1
   globeAtt@title = "MM5: daily averages"
   globeAtt@acronym = "MM5: Mesoscale Model"
   globeAtt@description = "Data analyzed uses archived NCEP reanalysis.
This project is in association with the NARCCAP project"
   globeAtt@process = "NetCDF is sorted output (top 10% extreme
events)"

   fileattdef(ncdf ,globeAtt)

end

-- 
Sho Kawazoe

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 12 16:06:43 2010