Re: convert spectral coefficients to grid

>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