Re: Problem of area weighted sum using "wgt_areasum2" function

For quick an dirty.
I extracted the Example 2 code, which assumes constant
lat spacing [it is not] ...

I calculated: 3.269624e+14. This is 91% of the wiki value.

If you account for variable lat spacing, this will increase.

I would look carefully at the masked grid points.

D

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************
begin
;************************************************
; Read in shferature
;************************************************
  ds = addfile("/project/cas/shea/uWisc_jun/ocn_g.nc","r")

  ar = ds->og
  lat = ds->lat
  lon = ds->lon
;*************************************************
; Area average
;*************************************************
    lat1 = new(dimsizes(lat({-90:90})),float)
    lon1 = new(dimsizes(lon({0:360})),float)
    lat1 = lat({-90:90})
    lon1 = lon({0:360})

    jlat = dimsizes( lat1 )
    rad = 4.0*atan(1.0)/180.0
    re = 6371220.0
    rr = re*rad
    dlon = abs(lon1(2)-lon1(1))*rr

    dx = dlon*cos(lat1*rad)
    dy = new ( jlat, typeof(dx))
    dy(0) = abs(lat1(2)-lat1(1))*rr
    dy(1:jlat-2) = abs(lat1(2:jlat-1)-lat1(1:jlat-2))*rr*0.5
    dy(jlat-1) = abs(lat1(jlat-1)-lat1(jlat-2))*rr

    area0 = dx*dy ; cell area function of latitude only
    area = new((/dimsizes(lat1),dimsizes(lon1)/),float)
    area = conform(area,area0,0)

    ari = wgt_areasum2(ar({-90:90},{0:360}), area, 0)

  print(ari)
; ===============================

   re = 6.37122e06
   rad = 4.0 * atan(1.0) / 180.
   con = re * rad
   clat = cos(lat * rad) ; cosine of latitude

   dlon = (lon(2) - lon(1)) ; assume dlon is constant
   dlat = (lat(2) - lat(1)) ; assume dlat is constant

   dx = con * dlon * clat ; dx at each latitude
   dy = con * dlat ; dy is constant
   dydx = dy * dx ; dydx(nlat)

   wgt = conform(ar, dydx, 0)

   ocn_ncl = wgt_areasum2(ar, wgt, 0) ; => qSum(ntim, klev)

   print(ocn_ncl)

   ocn_wiki = 3.61*10^14

   pc = (ocn_ncl/ocn_wiki)*100

   print(pc)

end

jun wrote:
> Hello
>
> When I using function "wgt_areasum2" to calculate the total sum of
> heat flux for each ocean, the results look like some weird, so I use
> the following script to test this function of area weighted sum. The
> variable "ar" in ocn_g.nc file has been set to "1" in ocean grid and
> to missing value at land grid. So the sum using this function should
> be the total area of global ocean.
>
> The calculating result of global ocean area is 1.6348*10^14*m^2, but
> the actual value should be 3.61*10^14*m^2 (refer
> to: http://en.wikipedia.org/wiki/Ocean).
>
> So, could you give me some suggestion about using of this function?
>
> Thanks.
>
> Jun Cheng
>
>
>
> ; Area of global ocn
> ;************************************************
> 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"
> ;************************************************
> begin
> ;************************************************
> ; Read in shferature
> ;************************************************
> ds = addfile("/ptmp/jcheng2/ocn_g.nc","r")
>
> ar = ds->og
> lat = ds->lat
> lon = ds->lon
> ;*************************************************
> ; Area average
> ;*************************************************
> lat1 = new(dimsizes(lat({-90:90})),float)
> lon1 = new(dimsizes(lon({0:360})),float)
> lat1 = lat({-90:90})
> lon1 = lon({0:360})
>
> jlat = dimsizes( lat1 )
> rad = 4.0*atan(1.0)/180.0
> re = 6371220.0
> rr = re*rad
> dlon = abs(lon1(2)-lon1(1))*rr
>
> dx = dlon*cos(lat1*rad)
> dy = new ( jlat, typeof(dx))
> dy(0) = abs(lat1(2)-lat1(1))*rr
> dy(1:jlat-2) = abs(lat1(2:jlat-1)-lat1(1:jlat-2))*rr*0.5
> dy(jlat-1) = abs(lat1(jlat-1)-lat1(jlat-2))*rr
>
> area0 = dx*dy ; cell
> area function of latitude only
> area = new((/dimsizes(lat1),dimsizes(lon1)/),float)
> area = conform(area,area0,0)
>
> ari = wgt_areasum2(ar({-90:90},{0:360}), area, 0)
>
> print(ari)
>
> end
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> 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
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk