Help to find daily total

Sir,
1.Version 5.1.1
2.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 "/home/akila/ncl/lib/scripts/my_scripts.ncl"

begin

   nlat = 94
   mlon = 192
   lon = ispan(0,mlon-1,1)*1.875-180
; lat = ispan(0,nlat-1,1)*1.875-88.54195
  print(lon)
;printVarSummary(lon)
  glat = gaus(nlat/2)
  lat =glat(:,0)
  lat=lat(::-1)
print(lat)

   lat!0 = "lat"
   lat&lat = lat
   lat_at_long_name = "latitude"
   lat_at_units = "degrees_north"
   lon!0 = "lon"
   lon&lon = lon
   lon_at_long_name = "longitude"
   lon_at_units = "degrees_east"

  ; miscellaneous
   tunits = "hours since 1990-04-01 00:00:0.0"
   minute = 0
   second = 0.d0
   filAtt = True
   nline = inttochar(10)

   lonout = ispan(0,mlon-1,1)*1.875 -180
   lonout!0 = "lon"
   lonout&lon = lonout
   lonout@long_name = "Longitude"
   lonout@units = "degrees_east"

   lat@long_name = "Latitude"
   lat@units = "degrees_north"

   mday = (/0,31,28,31,30,31,30,31,31,30,31,30,31/)

   diri = "/xxxxx/"
   diro = "/xxxxx/"

   p1 = new((/nlat,mlon/),float)

   do iyr = 1990,1990
    if(mod(iyr,400).eq.0 .or. mod(iyr,4).eq.0 .and. mod(iyr,100).ne.0)then
     mday(2) = 29
    else
     mday(2) = 28
    end if

    do imo = 04,12
     fout = diro+"prec.flx."+sprinti("%.2d",iyr)+sprinti("%.2d",imo)+".dy.nc"
     print(fout)

     time = new(mday(imo),double,"No_FillValue")
     time!0 = "time"
     time@units = "hours since 1800-04-01 00:00:00"
     time@long_name = "Time"

     prec = new((/mday(imo),nlat,mlon/),float,"No_FillValue")
     prec@_FillValue = 1e+20
     prec!0 = "time"
     prec!1 = "lat"
     prec!2 = "lon"
     prec&lat = lat
     prec&lon = lonout
     prec@long_name = "Precipitation"
     prec@units = "mm"
; prec@units = "kg/m^2/s"

     prec = 0.0
     
     do idy = 01, mday(imo)
; do idy =09,31
      time(idy-1) = ut_inv_calendar(iyr,imo,idy,12,0,0,time@units,0)

      do ihr = 00, 00
       fin = diri+"flx."+"fd"+sprinti("%d",iyr)+sprinti("%.2d",imo)+sprinti("%.2d",idy)+sprinti("%.2d",ihr)+".grb"
; print(iyr+" "+imo+" "+idy+" "+ihr+" "+fin)

       f = addfile(fin,"r")
; printVarSummary(f)
       fvname=getfilevarnames(f)

      vname1="PRATE_98_SFC_ave"
      vname2="PRATE_98_SFC_ave3h"
      vname3="PRATE_98_SFC_10"

     if(any(fvname .eq. vname1))then
        p0 = f->$vname1$
     end if
     if(any(fvname .eq. vname2))then
        p0 = f->$vname2$
     end if
     if(any(fvname .eq. vname3))then
        p0 = f->$vname3$
     end if

       p0!0 = "lat_98"
       p0!1 = "lon_98"
       p = p0(lat_98|:,lon_98|:)
       p1(:,0:mlon/2-1) = p(:,mlon/2:mlon-1)
       p1(:,mlon/2:mlon-1) = p(:,0:mlon/2-1)
       p1@_FillValue = -1e+20
       prec(idy-1,:,:) = prec(idy-1,:,:) + p1*3600*24

       delete(p0)
       delete(p)
       delete(f)
      end do
; id =id+1
     end do
     time&time = time
;print(time)
     prec&time = time
     printVarSummary(prec)
     printMinMax(prec,0)

     fnc = addfile(fout,"c")
     setfileoption(fnc,"DefineMode",True)

     fatt = True
     fatt@title = "grb flux data"
     fatt@Accumulation_Hour = "Daily Average"
     fileattdef(fnc,fatt)

     dimname = (/"time","lat","lon"/)
     dimsize = (/mday(imo), nlat, mlon/)
     dimunlim = (/True, False, False/)
     filedimdef(fnc,dimname,dimsize,dimunlim)

     filevardef(fnc, "time", typeof(time), getvardims(time))
     filevardef(fnc, "lat", typeof(lat), getvardims(lat))
     filevardef(fnc, "lon", typeof(lonout), getvardims(lonout))
     filevardef(fnc, "prec", typeof(prec), getvardims(prec))

     filevarattdef(fnc,"time",time)
     filevarattdef(fnc,"lat",lat)
     filevarattdef(fnc,"lon",lonout)
     filevarattdef(fnc,"prec",prec)

     setfileoption(fnc,"DefineMode",False)

     fnc->time = (/time/)
     fnc->lat = (/lat/)
     fnc->lon = (/lonout/)
     fnc->prec = (/prec/)

     delete(prec)
     delete(time)
   end do
   end do
end
3.want to calculate daily average
regards,
amm

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Dec 23 02:52:07 2009