>I'm interested in converting from spectral coefficients (T106) to a
>regular lat/lon grid. I assume NCL includes a command for this, since
>it includes Spherepack, but have not been able to find anything in the
>documentation. Any help would be great.
________________________
NCL can
[1] read the spectral coefficients.
[2] synthesize the coefficients to a gaussian grid
[3] interpolate to another grid.
I assume the original data are from a GRIB file?
If so, on the comman line use ncl_filedump
to examine the file's contents.
http://www.ncl.ucar.edu/Document/Tools/ncl_filedump.shtml
%> ncl_filedump foo.grb
Pick the variable u want to read/synthesize/regrid
load "./cc2g_time.ncl"
f = addfile ("foo.grb" , "r")
xg = cc2g (f, "TMP_GDS50_HYBL", (/128,256/), 0)
--- The cc2g_time is more complex than it needs to be. NCL can do what you want in a couple of steps. However, some grib files are huge and much memory is required. Hence, I broke it up such that it processed one time step and level at a time. The code returns a gaussian grid. You can [a] return a "regular" grid by commenting the "g2gsh" functions and replacing with [say] var(nt,:,:) = g2fsh(temp, (/nlat,mlon/)) or [b] xg = cc2g (f, "TMP_GDS50_HYBL", (/128,256/), 0) xf = g2fsh_Wrap(xg , (/181,360/)) good luck D g2gsh(temp, (/nlat,mlon/), tWave)
; --------------------------------------------------------
; Complex Coefficients to user defined Gaussian Grid
;
; The ascalar variable returned will either be a
; 3D (time,lat,lon) or 4D (time,lev,lat,lon) variable
; All time steps are done for the input variable.
; It is assumed that the latitudes are in ascending order.
;
; sample usage: x = cc2g (f, "TMP_GDS50_HYBL", (/64,128/), 42)
; sample usage: x = cc2g (f, vName, (/jlat,ilon/), tWave)
function cc2g (f:file, varName:string, dimOut[*], tWave:integer)
begin
temp = new ( (/160,320/), float) ; full array
a = new ( (/160,160/), float) ; real coef
b = new ( (/160,160/), float) ; imag coef
temp = 0.0 ; initilize
a = 0.0
b = 0.0
dimSizes = filevardimsizes(f,varName) ; get dimension sizes
numDims = dimsizes(dimSizes) ; number of dimensions
ntim = dimSizes(0) ; number of time steps
nlat = dimOut(0) ; # output latitudes
mlon = dimOut(1) ; # output longitudes
;print ("ntim="+ntim+" nlat="+nlat+" mlon="+mlon)
if (numDims.eq.5) ; (time,lev,2,lat,lon)
klvl = dimSizes(1) ; number of levels
;print (klvl)
var = new( (/ntim,klvl,nlat,mlon/), float) ; allocate space
do nt=0,ntim-1
do kl=0,klvl-1
;print ("nt="+nt+" kl="+kl)
a(0:106,0:106) = f->$varName$(nt,kl,0,:,:)
b(0:106,0:106) = f->$varName$(nt,kl,1,:,:)
shsgc(a,b,temp)
if (nlat.eq.160 .and. tWave.eq.0) then
var(nt,kl,:,:) = temp ; original grid
else
var(nt,kl,:,:) = g2gsh(temp, (/nlat,mlon/), tWave)
end if
end do
end do
else ; (time,lat,lon)
var = new( (/ntim,nlat,mlon/), float) ; allocate space
do nt=0,ntim-1
a(0:106,0:106) = f->$varName$(nt,0,:,:)
b(0:106,0:106) = f->$varName$(nt,1,:,:)
shsgc(a,b,temp)
if (nlat.eq.160 .and. tWave.eq.0) then
var(nt,:,:) = temp ; original grid
else
var(nt,:,:) = g2gsh(temp, (/nlat,mlon/), tWave)
end if
end do
end if
return(var)
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jun 21 2006 - 14:44:40 MDT
This archive was generated by hypermail 2.2.0 : Wed Jun 21 2006 - 17:13:25 MDT