computing daily precip from hourly data

Hi,

 I have a script that was discussed in a earlier talk thread (and is also
attached to this thread). I am looking for a way to adjust the starting time
with this script. Originally, this script reads in the entire file, and
converts hourly to daily from the first time step in that file. I am looking
at a time frame from 6UTC - 6UTC, so this script will need a offset.

 Comparing the raw input data with unit conversion (simply multiplied each 3
hr data by 10800, and summed up 8 time steps to get mm/day) with the script
using the entire file gives the same values However, trying to change the
ntStrt value, or doing a simple kitchen sink method to the input data (ie
start the file one time step further) to accommodate the offset gives
incorrect results.

Checking values was done with a simple kitchen sink method,
ncks -d time,0,40 -d xc,1,1 -d yc,1,1 .... to check values at a
certain location, for simple comparison purposes.

Any help is appreciated. I apologize if my question is a bit vague.

;*******************************************************************
; Load NCL library files.
;*******************************************************************

  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

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

;**********************************************************
; MAKE SURE you adjust the time the file reads (currently
; set at 3UTC - 3UTC) to time you want to analyze.
;**********************************************************

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

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

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

    ntJump = 8

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

    ntStrt = 0
    ntLast = ntJump-1
    do nt=0,NTIM-1,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)

    x@info_tag = "daily average"

;**************************************************************
; Convert units kg/(m2-s) to mm/day
; multiply by (10^3mm m-1 x 86400 s day-1) and
; divide by density_H2O (1000 kg m-3):
; [kg/m2-s][1000 mm/m][86400 s/day][(1/1000) m3/kg] ==> mm/day
;**************************************************************

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

;******************************************************************
; Create a new nc output file
; Also output variables directly into created netcdf file
;******************************************************************

   system("/bin/rm -f daily_output.nc") ; remove any pre-existing
file
   fout = addfile("daily_output.nc","c")
   filedimdef(fout,"time",-1,True) ; Makes time an UNLIMITED
dimension (as suggested by NCL)

   fout->pr = x
   fout->lat = lat
   fout->lon = lon
   fout->Lambert_Conformal = Lambert_Conformal

end

-- 
Sho Kawazoe

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 25 12:49:24 2010