writing netcdf

From: Ben Foster (foster AT hao.ucar.edu)
Date: Sun Oct 02 2005 - 14:25:31 MDT

Hi ncl:

I am reading a netcdf file containing 4-d model data (time,lev,lat,lon),
and writing a new one, equivalent to the input file except with the bottom
boundary extended downward several levels (6 new levels in this case).

My script (see attached) seems to be working ok so far, except that the
values on the output file of one of the two vertical coordinate variables
is incorrect. When I print the values of this coordinate var (lev(lev)) to
stdout from the script immediately after writing them to the output file,
they look fine, but an ncdump of the output file after script execution shows
wrong values. The other vertical coord var (ilev(ilev)) is correct on the
output file.

To find where the problem is occurring in the script, please search on
"Ask ncl" in the attached script.

My main problem is stated above, but as you can see, this is kind of a
brute force script, so I appreciate any suggestions about how to generalize
this, or make it more efficient.

BTW, in case you want to try running the script, I have put the netcdf
input file on http://download.hao.ucar.edu/pub/foster/for_ncl
(a copy of the script is there too).

Thanks for looking at this,


Ben Foster High Altitude Observatory (HAO)
foster@ucar.edu phone: 303-497-1595 fax: 303-497-1589
Nat. Center for Atmos. Res. P.O. Box 3000 Boulder CO 80307 USA

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/shea_util.ncl"
  infile_name = "FOSTER.tiegcm.ptest001.nc"
  ncin = addfile(infile_name ,"r")

  outfile_name = "FOSTER.tiegcm.newlbc.nc"
  system("rm "+outfile_name)
  ncout = addfile(outfile_name,"c")

  varnames = getfilevarnames(ncin)
  varnames = varnames(::-1) ; reverse order
  nvars = dimsizes(varnames)
  nlev = getfilevardimsizes(ncin,"lev") ; number of levels on input file
; print("nvars="+nvars+" varnames="+varnames+" nlev="+nlev)

  newnlev = 35
; number of new lower levels (old lbc downward to new lbc)
; (assumes newnlev > nlev)
  nnewlev = newnlev-nlev
; New interface array:
  newilev = new(newnlev,"double")
  newilev = (/ -10., -9.5, -9., -8.5, -8., -7.5, -7., -6.5, -6., -5.5, -5., -4.5, -4., -3.5, -3., -2.5, -2., -1.5, -1., -0.5, 0., 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5, 7. /)
; New midpoints array (half dz above interfaces, assumes dz of 0.5):
  newlev = newilev+0.25
; Define new vertical dimension sizes:
  filedimdef(ncout,"lev" ,newnlev,False)

; print("newilev="+newilev)
; print("newlev="+newlev)
; Variable loop:
; do i=0,nvars-1 ; usually around 67 for tiegcm primary history file
  do i=0,40 ; catches a few 4d vars
    f = ncin->$varnames(i)$
    vdims = getfilevardims(ncin,varnames(i))
    vdimsizes = getfilevardimsizes(ncin,varnames(i))
; vattnames = getfilevaratts(ncin,varnames(i)) ; same as getvaratts??
    vatts = getvaratts(f) ; for info only
; print("i="+i+" var "+varnames(i)+" vdims="+vdims+" vdimsizes="+vdimsizes)
; Make new var if current var has a lev dimension:
; (Current var should have either 1 (coord vars lev and ilev), or 4 dimensions)
; Might want to ask ncl people how to generalize this.
    if (any(vdims.eq."lev").or.any(vdims.eq."ilev")) then
; print(" ")
; printVarSummary(f)
      vtype = getfilevartypes(ncin,varnames(i))
      ndims = dimsizes(vdims)
; Make new 1-d coord var (either lev(lev) or ilev(ilev)): <-- need to add mlev
      if (ndims.eq.1) then ; 1-d coord var (either lev(lev) or ilev(ilev)
        print("Make new 1-d var "+varnames(i)+" vdims="+vdims+" vtype="+vtype)
        print(" vatts = "+vatts)
        if (varnames(i).eq."lev") then ; midpoints coord var
          fnew = newlev ; is a double
        else ; interface coord var
          fnew = newilev ; is a double
        end if
        print("assign values for new var "+varnames(i)+": vdims="+vdims+" fnew="+fnew)
; Make new 4-d coord var (time, lev, lat, lon):
      else ; this "else" makes the big assumption that this is a 4d field
        if (ndims.eq.4) then
          print("Make new 4-d var "+varnames(i))
          print(" vatts = "+vatts)
          fnew = new((/vdimsizes(0),newnlev,vdimsizes(2),vdimsizes(3)/),vtype)
; Assign fnew from f for old levs, leaving missing values from new lbc up to old lbc:
          fnew(:,nnewlev:newnlev-1,:,:) = f(:,0:nlev-1,:,:)
; Assign values for current 4d var at new levels here. For now, use the old lbc:
; (could put zeroes here, or leave it _FillValue, or even import/read data)
; fnew(:,0:nnewlev-1,:,:) = f(:,0,:,:) ; ncl will not do this
          do k=0,nnewlev-1 ; this does not slow it down too much
            fnew(:,k,:,:) = f(:,0,:,:)
          end do
          print(">>> Did not expect to find ndim="+ndims+" for var "+varnames(i))
        end if
      end if ; 1d or 4d var
; Write new var to file:
      ncout->$varnames(i)$ = fnew

; Ask ncl: This prints to stdout ok, but lev(lev) values are wrong on output file!
      if (varnames(i).eq."lev") then
        print("Wrote new levs: fnew="+fnew)
      end if

; Clean out ncl var fnew for next var:
; Copy attributes (assuming no changes for now, 10/2/05)
; Give new value to LBC:
; (This else is a big assumption)
    else ; var does not have a lev dimension
      if (varnames(i).eq."LBC") then
        f = (/ newilev(0) /)
      end if
      ncout->$varnames(i)$ = f
    end if
; Clean up current var:
  end do ; i=0,nvars-1


ncl-talk mailing list

This archive was generated by hypermail 2b29 : Mon Oct 03 2005 - 11:00:18 MDT