Re: Problem computing grid points

Brandon Fisel wrote:
> Dear NCL'ers
>
> I'm attempting to compute two sums below, poleward of a specified
> latitude and plot them using a time series. The .nc file contains
> daily averaged sea ice fractions for 1998.
> Sum 1 = ( # of grid points with 50% < fraction < 90% ) x cos(lat)
> Sum 2 = ( all grid points x cos(lat))
The code you provided is a bit convoluted. I had trouble following it.
Try the following:

  fsiLo = 50 ; %
  fsiHi = 90 ; %
  latLo = 70
  latHi = 90

  diri = "./"
  fili = "met_em.1998.seaice.nc"
;*****************************************
; Data input
;*****************************************
  f = addfile(diri+fili, "r")
  fsi = f->fsi ; ( Time, south_north, west_east )
  xlat = f->XLAT_M ; ( south_north, west_east )
  Times_chr= f->Times
  Times_str= chartostring( Times_chr)
  ntim = dimsizes(Times_str)

;*****************************************
; Compute weighted averages
;*****************************************
  clat = where(xlat.ge.latLo .and. xlat.le.latHi \
                  ,cos(0.01745329*xlat), 1e20)
  clat@_FillValue = 1e20
  clat_sum = sum(clat)

  wfsi = new ( ntim, typeof(fsi), getFillValue(fsi))
  do nt=0,ntim-1
     wfsi(nt) = sum(fsi(nt,:,:)*clat)/clat_sum
  end do
  wfsi_at_long_name = "seaice concentration"
  wfsi_at_units = "%"

>
> Below is part of the script:
>
> ; Add file and grab attributes
> f = addfile("met_em.1998.seaice.nc <http://met_em.1998.seaice.nc>","r")
> print(f)
> ; Get times for the file and set attributes
> TimeChar = f->Times
> DimTimeChar = dimsizes(TimeChar)
> nTime = DimTimeChar(0)-1 ; convenience
> times = wrf_user_list_times(f)
> ntimes = dimsizes(times)
> print ("wrf_user_list_times is:" + ntimes)
> lat = f->XLAT_M(:,:)
> DimLat = dimsizes(lat)
> nS_N = DimLat(0)
> nW_E = DimLat(1)
> lat_at_long_name = "Latitude"
> lat_at_standard_name = "latitude"
> lat_at_units = "degrees latitude"
> lat!0 = "south_north"
> lat!1 = "west_east"
> lon = f->XLONG_M(:,:)
> lon_at_long_name = "Longitude"
> lon_at_standard_name = "longitude"
> lon_at_units = "degrees longitude"
> lon!0 = "south_north"
> lon!1 = "west_east"
> ; Plot daily critical fractions
> do nt=0,nTime
> type = "ps"
> prefix = "daily_cri_fra_"
> suffix = times(nt)+""
> print ("Suffix is:" + suffix)
> path = diro+""+prefix+""+suffix
> print ("Path is:" + path)
> wks = gsn_open_wks(type,path)
> dims = getfilevardims(f,"fsi")
> si = f->fsi(nt,:,:)
> lat2d = lat
> lon2d = lon
> si_at_lat2d = lat2d
> si_at_lon2d = lon2d
> si = si/100
> printVarSummary(si)
> si_min = .50 ; minimum fraction
> si_max = .90 ; maximum fraction
> ; Not certain from here on
> lats =
> new((/75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90/),float) ;
> create array of lats
> do nLat=75,lats-1
> ice(nLat) = si(lat|:,lon|:)
> cr_level = dim_num(ice.ge.si_min .and. ice.le.si_max) ; type
> integer sum of grid points
> copy_VarCoords(ice,cr_level)
> sum_crlv = cr_level*cos(nLat)
> sum_tot = sum_crlv(nLat)
> end do
>
> Thank you in advance!
>
> -------------------------------------------------------------
> Brandon Fisel
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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