Some useful functions:
http://www.ncl.ucar.edu/Document/Functions/Contributed/lonGlobeF.shtml
http://www.ncl.ucar.edu/Document/Functions/Contributed/latGau.shtml
===
Curious .... Why do you have
mday = (/0,31,28,31,30,31,30,31,31,30,31,30,31/)
Why is the zero present?
===
http://www.ncl.ucar.edu/Document/Functions/Built-in/days_in_month.shtml
is better/cleaner than manual specification. In your case
mday = days_in_month( iyr, imo )
====
Why?
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)
===
do idy = 01, mday(imo)
time(idy-1) = ut_inv_calendar(iyr,imo,idy,12,0,0,time@units,0)
do ihr = 00, 00
This is only one hour??? Why 00? I assume it should be something like
do ihr = 0, 21,3
or
do ihr = 0, 18,6
*****************************************************************
Recommendations:
[0] I would rewrite the code :-)
[1] For an example, use the following and compare outputs
%> ncl_filedump flx.....grb | less
There will be no "time" dimension. Next, try
%> ncl_filedump -itime flx.....grb | less
You will see a "time" dimension ["initial_time0_hours"]
float FOO ( initial_time0_hours, lat_98, lon_98 )
You will also see that NCL has provided
(a) the latitudes and longitudes. *No* need for you to compute them.
(b) the "time" variable ... *No* need for you to compute them.
Use
http://www.ncl.ucar.edu/Document/Functions/Built-in/setfileoption.shtml
Specifically,
setfileoption("grib", "SingleElementDimensions", "All")
to make NCL include a time dimension. Place this before reading any
grib files.
[2] Use addfiles
http://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml
Then I would read all the files for a specific day.
fils = systemfunc("cd "+diri+" ; flx.fd"+
sprinti("%d",iyr)+sprinti("%.2d",imo)+sprinti("%.2d",idy)+"*")
print(fils)
f = addfiles(diri+fils+.grb", "r")
p0 = f[:]->$vname$ ; read all time steps for this
day (time,lat,lon)
psum = dim_sum_n_Wrap(p0, 0) ; sum all the hourly values
psum!0 = "time"
psum!1 = "lat"
psum!2 = "lon"
psum@long_name = "daily Precip Total"
psum@units = "mm"
printVarSummary(psum)
[4] Use the simple method to write the netCDF. For small
files this is easier/cleaner. Method 1
http://www.ncl.ucar.edu/Applications/o-netcdf.shtml
Good luck
===============
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
lat = latGau(nlat, "lat", "latitude", "degrees_north")
printVarSummary(lat)
lon = lonGlobeF(mlon, "lon", "longitude", "degrees_east")
printVarSummary(lon)
; miscellaneous
tunits = "hours since 1990-04-01 00:00:0.0"
minute = 0
second = 0.d0
filAtt = True
nline = inttochar(10)
diri = "/xxxxx/"
diro = "/xxxxx/"
do iyr = 1990,1990
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 = "Tim
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
On 12/23/2009 02:52 AM, Ammu Priya wrote:
> 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
>
-- ====================================================== Dennis J. Shea tel: 303-497-1361 | P.O. Box 3000 fax: 303-497-1333 | Climate Analysis Section | Climate& Global Dynamics Div. | National Center for Atmospheric Research | Boulder, CO 80307 | USA email: shea 'at' ucar.edu | ======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Dec 23 14:23:56 2009
This archive was generated by hypermail 2.1.8 : Tue Dec 29 2009 - 10:29:16 MST