# Re: Help to find daily total

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Wed Dec 23 2009 - 14:22:39 MST

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?
===

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.

Specifically,
setfileoption("grib", "SingleElementDimensions", "All")
to make NCL include a time dimension. Place this before reading any
grib files.

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)

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

Good luck

===============

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)

; 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)

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
>
> 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)
>
> ; 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)
>
> 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