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-talkReceived on Thu Apr 30 2009 - 10:21:41 MDT
This archive was generated by hypermail 2.2.0 : Thu Apr 30 2009 - 15:56:01 MDT