Fw: convert grb to nc

----- Forwarded Message ----
From: Ammu Priya <ammu9016@ymail.com>
To: ncl-talk@ucar.edu
Sent: Wed, December 23, 2009 3:22:01 PM
Subject: 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 convert grb giles to nc files
4.example filesflx.fd1990072600.grb flx.fd1990121312.grb sigp.fd1990102000
flx.fd1990072603.grb
flx.fd1990072606.grb flx.fd1990121318.grb sigp.fd1990111900
flx.fd1990072609.grb flx.fd1990121321.grb sigp.fd1990120400
flx.fd1990072612.grb flx.fd1990121400.grb sigp.fd1990121900
flx.fd1990072615.grb flx.fd1990121403.grb sigp.fd1991010300
flx.fd1990072618.grb flx.fd1990121406.grb sigp.fd1991011800
flx.fd1990072621.grb

regards,
amm

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