Re: write the file into netcdf

Hi Mary,

It's a bit hard to debug your code because you have several "do" loops and they are indented somewhat randomly. It's hard to tell which code is inside which do loop.

I think your problem is with this code:

> do x = 0,405
> do y = 0,415
> . . .
> end do
> x = x +1
> end do
> y = y +1
> end do

You shouldn't increment x and y by 1 before the "end do", because the "do x=0,405 and "do y=0,415" is doing this for you.

Remove both of these lines:

> x = x +1
> y = y +1

There are some other things you can clean up, too. "do" loops can be very expensive, so try to keep redundant code *outside* of a do loop.

For example, you have this code inside the do loop:

> Tmin= new((/31,531,420/), float)
> Tmax = new((/31,531,420/), float)
> rain = new((/31,531,420/), float)
> sw = new((/31,531,420/), float)

You can move this outside the loop, so these variables are only created once.

Also, for getting the file names, do you really need all this code?

> do i = 1,12
> if (i .lt. 10) then
> tmp_files = systemfunc("ls /2001/20010" + i + "*LDASIN*.nc")
> else
> tmp_files = systemfunc("ls /2001/2001" + i + "*LDASIN*.nc")
> end if
> files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> n1 = n1 + dimsizes(tmp_files)
> print(n1)
> delete(tmp_files)
> end do

I would think you could just do this:

     files = systemfunc("ls /2001/2001[0-1]*LDASIN*.nc")

Finally, it's a good idea to not not hard-code dimension sizes, if possible.

So, instead of this:

> do x = 0,405
> do y = 0,415

you can do this:

;---Get dimension sizes.
dims = dimsizes(f[:]->T2D)
nv2 = dims(0)
nx = dims(1)
ny = dims(2)
do x = 0,nx-1
do y = 0,ny-1

--Mary

On Nov 18, 2012, at 1:13 PM, mary charusombat wrote:

> Dear all (I do not know I should reply to only Mary or all)
>
> From my script, I changed the switch x and y value (as Mary recommend) but it still does not work.
> Then I try a coupled grid number and found out that I can not go beyond 405 in x direction and 415 in y direction.
> The original image size is 531 x 419.
> Is it ncl bug or the limitation in ncl?
>
> And when I wrote to new nc files.
> It has no values in there.
> I got the header of my new file like this.
> How I write into new netcdf with lat and lon, my parameter and time (it should be 31 (ncl1)?
>
> netcdf \2001_daily {
> dimensions:
> time = UNLIMITED ; // (0 currently)
> ncl1 = 31 ;
> ncl2 = 531 ;
> ncl3 = 420 ;
> ncl4 = 31 ;
> ncl5 = 531 ;
> ncl6 = 420 ;
> ncl7 = 31 ;
> ncl8 = 531 ;
> ncl9 = 420 ;
> ncl10 = 31 ;
> ncl11 = 531 ;
> ncl12 = 420 ;
> variables:
> float Tmin(ncl1, ncl2, ncl3) ;
> Tmin:_FillValue = -999.f ;
> float Tmax(ncl4, ncl5, ncl6) ;
> Tmax:_FillValue = -999.f ;
> float sw(ncl7, ncl8, ncl9) ;
> sw:_FillValue = -999.f ;
> float rain(ncl10, ncl11, ncl12) ;
> rain:_FillValue = -999.f ;
>
>
>
> ________________________________________________
>
> ;*****************************************************
> ;xy_18.ncl
> ;*****************************************************
> 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 "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> ;*****************************************************
> begin
> delimiter = " "
> ;files = new(744, string)
> files = new(8760, string)
> n1 = 0
> a = "geo_em.d01.nc"
> do i = 1,12
> if (i .lt. 10) then
> tmp_files = systemfunc("ls /2001/20010" + i + "*LDASIN*.nc")
> else
> tmp_files = systemfunc("ls /2001/2001" + i + "*LDASIN*.nc")
> end if
> files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> n1 = n1 + dimsizes(tmp_files)
> print(n1)
> delete(tmp_files)
> end do
> f = addfiles(files, "r")
> a1 = addfile(a, "r")
> do x = 0,405
> do y = 0,415
> v2tem = f[:]->T2D(:,x,y)
> v3tem = f[:]->SWDOWN(:,x,y)
> v4tem = f[:]->RAINRATE(:,x,y)
> nv2 = dimsizes(v2tem)
> print(nv2)
> v2tem@_FillValue = -1e+36
> v3tem@_FillValue = -1e+36
> v4tem@_FillValue = -1e+36
> im = 0 ;start at 0700 UTC for MST time 24*14 = 336
> Tmin= new((/31,531,420/), float)
> Tmax = new((/31,531,420/), float)
> rain = new((/31,531,420/), float)
> sw = new((/31,531,420/), float)
> do i = 0,30
> Tmin(i,x,y) = min(v2tem(im:im+23))
> Tmax(i,x,y) = max(v2tem(im:im+23))
> sw(i,x,y) = sum(v3tem(im:im+23))
> rain(i,x,y) = sum(v4tem(im:im+23))
> im = im+24
> end do
> x = x +1
> end do
> y = y +1
> end do
> system ("/bin/rm -f 2001_daily.nc")
> ncdf = addfile("2001_daily.nc","c")
> filedimdef(ncdf,"time",-1,True)
> ncdf->Tmin = Tmin
> ncdf->Tmax = Tmax
> ncdf->sw = sw
> ncdf->rain = rain
> end
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> On Fri, Nov 9, 2012 at 1:15 PM, mary charusombat <marycharusombat@gmail.com> wrote:
> Dear all
> I try to do daily average of the hourly netcdf files.
> I wrote the scripts to extract data in each point and then average it in 24 hr for every points.
> Then I write the daily data back into a new netcdf file.
> I am not sure it is a good way to do that.
> This is my code and I get the script out of range and memory error.
> My netcdf files have
>
> // global attributes:
> :TITLE = "OUTPUT FROM CONSOLIDATE_GRIB v20110427" ;
> :missing_value = -1.e+36f ;
> :_FillValue = -1.e+36f ;
> :WEST-EAST_GRID_DIMENSION = 531 ;
> :SOUTH-NORTH_GRID_DIMENSION = 420 ;
> :DX = 4000.f ;
> :DY = -999.9f ;
> :TRUELAT1 = 42.627f ;
> :TRUELAT2 = 42.627f ;
> :LA1 = 34.48816f ;
> :LO1 = -103.878f ;
> :STAND_LON = -92.409f ;
> :MAP_PROJ = 1 ;
> :MMINLU = "USGS"
>
>
> Please help
> Mary
>
>
>
>
> This is my code
>
> begin
> delimiter = " "
> files = new(8760, string)
> n1 = 0
> a = "geo_em.d01.nc"
> do i = 1,12
> if (i .lt. 10) then
> tmp_files = systemfunc("ls ./2001/20010" + i + "*LDASIN*.nc")
> else
> tmp_files = systemfunc("ls ./2001/2001" + i + "*LDASIN*.nc")
> end if
> files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> n1 = n1 + dimsizes(tmp_files)
> print(n1)
> delete(tmp_files)
> end do
> f = addfiles(files, "r")
> a1 = addfile(a, "r")
> tmp_files = systemfunc("ls /scratch/lustreC/u/ucharuso/hrldas-v3.3/HRLDAS_COLLECT_DATA/run/2001/2001" + i + "*LDASIN*.nc")
> end if
> files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> n1 = n1 + dimsizes(tmp_files)
> print(n1)
> delete(tmp_files)
> end do
> f = addfiles(files, "r")
> a1 = addfile(a, "r")
> do x = 0,419
> do y = 0,530
> v2tem = f[:]->T2D(:,x,y)
> v3tem = f[:]->SWDOWN(:,x,y)
> v4tem = f[:]->RAINRATE(:,x,y)
> nv2 = dimsizes(v2tem)
> print(nv2)
> v2tem@_FillValue = -1e+36
> v3tem@_FillValue = -1e+36
> v4tem@_FillValue = -1e+36
> im = 0 ;start at 0700 UTC for MST time 24*14 = 336
> Tmin= new((/365,420,531/), float)
> Tmax = new((/365,420,531/), float)
> rain = new((/365,420,531/), float)
> sw = new((/365,420,531/), float)
> do i = 0,364
> Tmin(i,x,y) = min(v2tem(im:im+22))
> Tmax(i,x,y) = max(v2tem(im:im+22))
> sw(i,x,y) = sum(v3tem(im:im+22))
> rain(i,x,y) = sum(v4tem(im:im+22))
> im = im+22
> end do
> x = x +1
> end do
> y = y +1
> end do
> system ("/bin/rm -f 2001_daily.nc")
> ncdf = addfile("2001_daily.nc","c")
> filedimdef(ncdf,"time",-1,True)
> ncdf->Tmin = Tmin
> ncdf->Tmax = Tmax
> ncdf->sw = sw
> ncdf->rain = rain
> end
>
> fatal:Subscript out of range, error in subscript #1
> ^Mfatal:Memory allocation failure:[errno=12]
> ^Mfatal:Execute: Error occurred at or near line 30 in file gridextract.ncl
>
>
>
> On Fri, Nov 9, 2012 at 4:57 PM, Mary Haley <haley@ucar.edu> wrote:
> The error is telling you that a subscript is out of range.
>
> From the script you included below, line 30 would be this line:
>
> v2tem = f[:]->T2D(:,x,y)
>
> My guess is that either "x" or "y" are not in the range of the size of the "T2D" variable on the file.
>
> Do you perhaps have these swapped? That is, should it be:
>
> v2tem = f[:]->T2D(:,y,x)
>
> Are you getting any other errors out of this script? It looks like you have an "end if" statement that has no corresponding "if" statement at about line 19.
>
> --Mary
>
> On Nov 9, 2012, at 11:15 AM, mary charusombat wrote:
>
> > Dear all
> > I try to do daily average of the hourly netcdf files.
> > I wrote the scripts to extract data in each point and then average it in 24 hr for every points.
> > Then I write the daily data back into a new netcdf file.
> > I am not sure it is a good way to do that.
> > This is my code and I get the script out of range and memory error.
> > My netcdf files have
> >
> > // global attributes:
> > :TITLE = "OUTPUT FROM CONSOLIDATE_GRIB v20110427" ;
> > :missing_value = -1.e+36f ;
> > :_FillValue = -1.e+36f ;
> > :WEST-EAST_GRID_DIMENSION = 531 ;
> > :SOUTH-NORTH_GRID_DIMENSION = 420 ;
> > :DX = 4000.f ;
> > :DY = -999.9f ;
> > :TRUELAT1 = 42.627f ;
> > :TRUELAT2 = 42.627f ;
> > :LA1 = 34.48816f ;
> > :LO1 = -103.878f ;
> > :STAND_LON = -92.409f ;
> > :MAP_PROJ = 1 ;
> > :MMINLU = "USGS"
> >
> >
> > Please help
> > Mary
> >
> >
> >
> >
> > This is my code
> >
> > begin
> > delimiter = " "
> > files = new(8760, string)
> > n1 = 0
> > a = "geo_em.d01.nc"
> > do i = 1,12
> > if (i .lt. 10) then
> > tmp_files = systemfunc("ls ./2001/20010" + i + "*LDASIN*.nc")
> > else
> > tmp_files = systemfunc("ls ./2001/2001" + i + "*LDASIN*.nc")
> > end if
> > files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> > n1 = n1 + dimsizes(tmp_files)
> > print(n1)
> > delete(tmp_files)
> > end do
> > f = addfiles(files, "r")
> > a1 = addfile(a, "r")
> > tmp_files = systemfunc("ls /scratch/lustreC/u/ucharuso/hrldas-v3.3/HRLDAS_COLLECT_DATA/run/2001/2001" + i + "*LDASIN*.nc")
> > end if
> > files(n1: n1 + dimsizes(tmp_files) - 1) = tmp_files
> > n1 = n1 + dimsizes(tmp_files)
> > print(n1)
> > delete(tmp_files)
> > end do
> > f = addfiles(files, "r")
> > a1 = addfile(a, "r")
> > do x = 0,419
> > do y = 0,530
> > v2tem = f[:]->T2D(:,x,y)
> > v3tem = f[:]->SWDOWN(:,x,y)
> > v4tem = f[:]->RAINRATE(:,x,y)
> > nv2 = dimsizes(v2tem)
> > print(nv2)
> > v2tem@_FillValue = -1e+36
> > v3tem@_FillValue = -1e+36
> > v4tem@_FillValue = -1e+36
> > im = 0 ;start at 0700 UTC for MST time 24*14 = 336
> > Tmin= new((/365,420,531/), float)
> > Tmax = new((/365,420,531/), float)
> > rain = new((/365,420,531/), float)
> > sw = new((/365,420,531/), float)
> > do i = 0,364
> > Tmin(i,x,y) = min(v2tem(im:im+22))
> > Tmax(i,x,y) = max(v2tem(im:im+22))
> > sw(i,x,y) = sum(v3tem(im:im+22))
> > rain(i,x,y) = sum(v4tem(im:im+22))
> > im = im+22
> > end do
> > x = x +1
> > end do
> > y = y +1
> > end do
> > system ("/bin/rm -f 2001_daily.nc")
> > ncdf = addfile("2001_daily.nc","c")
> > filedimdef(ncdf,"time",-1,True)
> > ncdf->Tmin = Tmin
> > ncdf->Tmax = Tmax
> > ncdf->sw = sw
> > ncdf->rain = rain
> > end
> >
> > fatal:Subscript out of range, error in subscript #1
> > ^Mfatal:Memory allocation failure:[errno=12]
> > ^Mfatal:Execute: Error occurred at or near line 30 in file gridextract.ncl
> > _______________________________________________
> > ncl-talk mailing list
> > List instructions, subscriber options, unsubscribe:
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Nov 21 11:15:47 2012