Re: nco

From: Dennis Shea (shea AT XXXXXX)
Date: Thu Oct 26 2000 - 16:53:42 MDT

  • Next message: Jan Dutton: "4.2.2 and NCL error"

    Off line, I was sent the following:

    > My grid is global. It does not include pole points. The data contains
    > missing values.
    >
    > This is what I want to do.
    > 1. The data is at 0.5deg by 0.5deg resolution. (lat\lon).
    >
    > 2. I wish to compare it to another dataset which is at 2.5deg by 2.5deg
    > resolution.
    >
    > 3. Thus I want to reduce the resolution of data1 to 2.5 by 2.5deg
    >
    > 4. I want to do this by averaging over 5grid cells (lat) and 5 grid cells
    > (lon) so I get a value which will be representative of the area covered by
    > 2.5deg by 2.5deg grid cell.
    >
    > 5. I do not want to select every other 5th cell.
    >
    > 6. there are 20 files (20 years of data), each containing 12 months.

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

    I wrote two simple NCL functions which must be loaded prior to
    your main script. Pls note the use of NCL built-in functions
    and the array syntax. NCL automatically ignore missing values
    as indicated by _FillValue. I also include two test 'main' scripts.

    function CellAverage (x[*][*]:numeric, nbx:integer, nby:integer)
    ; compute "block" arithmetic averages
    ; missing values [ie: x@_FillValue ] automatically ignored
    ; Usage: xNew = CellAverage (x,nbx,nby)
    begin
      dimx = dimsizes(x) ; input dimension sizes
      jlat = dimx(0)
      ilon = dimx(1)

      nlat = jlat/nby ; output dimension sizes
      mlon = ilon/nbx
      y = new ( (/nlat,mlon/), typeof(x) )
                                ; pure arithmetic average
      nl = -1
      do jl=0,jlat-1,nby
         nl = nl+1
         ml = -1
        do il=0,ilon-1,nbx
           ml = ml+1
           y(nl,ml) = avg(x(jl:jl+nby-1,il:il+nbx-1)) ; NCL built-in func
        end do
      end do

      return (y)
    end
    ; ----------------------------------------------------------------
    function CellAverageCosWgt (x[*][*]:numeric, nbx:integer, nby:integer,
    lat[*]:numeric)
    ; compute "block" cosine wgted averages
    ; missing values [ie: x@_FillValue ] automatically ignored
    ; Usage: xNew = CellAverageCosWgt (x,nbx,nby)
    begin
      dimx = dimsizes(x) ; input dimension sizes
      jlat = dimx(0)
      ilon = dimx(1)

      rad = 4.0*atan(1.0)/180.
      clat = cos(lat*rad)

      wgt = new ( (/jlat,ilon/), float) ; use 2D array for conformity
      do il=0,ilon-1
         wgt(:,il) = clat
      end do

      nlat = jlat/nby ; output dimension sizes
      mlon = ilon/nbx
      y = new ( (/nlat,mlon/), typeof(x) )
                                ; cosine wgted average: sum(x*wgt)/sum(wgt)
                                ; sum is an NCL built-in function
      nl = -1
      do jl=0,jlat-1,nby
         nl = nl+1
         ml = -1
        do il=0,ilon-1,nbx
           ml = ml+1
           y(nl,ml) = sum(x(jl:jl+nby-1,il:il+nbx-1)*wgt(jl:jl+nby-1,il:il+nbx-1))/
    \
                      sum(wgt(jl:jl+nby-1,il:il+nbx-1))
        end do
      end do

      return (y)
    end

    ;-------------------- main 1
    begin
                         ; bogus input 0.5 degree grid
      jlat = 360
      ilon = 720
      x = new ( (/jlat,ilon/), float)
      x = 1000. ; set all value to 1000.

      nbx = 5 ; # of bins in lon [x] dir
      nby = 5 ; # of bins in lat [y] dir

      xNew = CellAverage (x, nbx, nby)

      print (" ============wgt by cos(lat)================== ")

      lat = ispan (0,jlat-1,1)*0.5 - 89.5 ; generate latitudes

      xNew = CellAverageCosWgt (x, nbx, nby, lat)
    end

    ;-------------------- main 2: related to (6) above
    begin
      nmos = 12
                         ; bogus input 0.5 degree grid
      jlat = 360
      ilon = 720
      x = new ( (/nmos,jlat,ilon/), float)
      x = 1000. ; set all value to 1000.

      nbx = 5 ; # of bins in lon [x] dir
      nby = 5 ; # of bins in lat [y] dir

      nlat = 72 ; jlat/nby
      mlon = 144 ; ilon/nbx
      xNew = new ( (/nmos,nlat,ilon/), typeof(x)) ; preallocate
      
      do nmo=0,nmos-1
         xNew(nmo,:,:) = CellAverage (x(nmo,:,:), nbx, nby)
      end do
    end



    This archive was generated by hypermail 2b29 : Fri Oct 27 2000 - 07:42:15 MDT