Visualizing data in Cubed-Sphere grid

IN PROGRESS

Cubed-Sphere grid background

The GEOS-5 now has products that are stored natively in cubed-sphere grid. The utilities listed here may help users visualize the data in cubed-sphere grid or map the data to lat-lon grid.

Fortran

To map the cubed-sphere grid data to lat-lon gird data, the program c2l_CFIO_offline.x is built with AGCM at src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/FVdycoreCubed_GridComp. The example configure file c2l_CFIO_offline.rc is also generated. It should be run as mpirun –np 6 ./c2l_CFIO_offline.x. This will produce a lat-lon grid data file whose name is configured in c2l_CFIO_offline.rc. Users then can use their familiar tools to visualize the lat-lon grid data.

Matlab

Users can find a “GridUtils” package in MATLAB ( R2015b and R2016a) at http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_HP_Output_Data#MATLAB_.28Cubed-sphere_and_regular_data.29 . The instructions are straight forward. Although it is intended for GCHP, it can be used for all the other outputs that use the same formats as that of GCHP.

IDL

The script read_and_interpolate_cube.pro can read data in the cubed-sphere grid in NetCDF files, interpolate the data , map the data in lat-lon format and plot. The script cube_to_latlon.pro uses read_and_interpolate_cube.pro to convert the variables in cubed-sphere grid to lat-lon grid and write the variables to NetCDF files. It should be noted that the tile files, input and output files are hard-wired in the code. The users will need to provide those files and change the code accordingly.

; read_and_interpolate_cube.pro
pro read_and_interpolate_cube, myfile, myVar, data, gLons, gLats, LonBeg, LonEnd, LatBeg, LatEnd, undef, TRIANGLES=TR, SCALE_FACTOR=SCALE_FACTOR, PIXEL=PIXEL, FLIP=FLIP, regridCubeToLatlon=regridCubeToLatlon, TIME=TIME, LEVEL=LEVEL, DATA_LEVELS=DATA_LEVELS, RANGE=RANGE, regridResolution=regridResolution, LOG=LOG, NCOLORS=NCOLORS, PROJECT=PROJECT, LISTVARS=LISTVARS, CMORPH_DATA=CMORPH_DATA, TRMM_DATA=TRMM_DATA, FILEID=INFILEID, STRETCH=STRETCH, CONSERV=CONSERV

 print, myfile
 help, myVar
;help, LonBeg, LonEnd, LatBeg, LatEnd
if (~keyword_set(TIME)) then TIME=0

exists = FILE_TEST(myfile)
if (exists) then begin
  if ( keyword_set(INFILEID)) then fileID = INFILEID
  if (~keyword_set(INFILEID)) then fileID = ncdf_open(myfile, /NOWRITE)

  if (keyword_set(LISTVARS)) then begin
    ;find the number of variables
    fileinq_struct = ncdf_inquire(fileID)
    nvars = fileinq_struct.nvars
    ; print output so far
    print, myfile
    print, nvars, ' variables in file'
    print, ' '
    print, ' '
    ;loop through variables
    for varndx = 0,nvars-1 do begin
     ; get the name, datatype, dims, number of attributes
       varstruct = ncdf_varinq(fileID,varndx)
     ; print the variable index, name, datatype, dims
       print, '--------------------------------------------------------'
       print,varndx,'    ',varstruct.name, '     ',varstruct.datatype
       print,'dims = ',varstruct.dim
    endfor ; variable loop
    stop
  endif

  finq = NCDF_INQUIRE( fileID )         ; finq /str = ndims, nvars, ngatts, recdim
  dim_unlimited = finq.recdim           ; = -1 if undefined, otherwise index into dim array
  if ( finq.ndims ge 2 ) then begin
   dimstr = ' '
   dimsize = 0L
   NCDF_DIMINQ, fileID, 0, dimstr, dimsize
  ;help, dimstr
   NX=dimsize
   NCDF_DIMINQ, fileID, 1, dimstr, dimsize
  ;help, dimstr
   NY=dimsize
  endif
  NZ=0L
  NT=1L
  if ( finq.ndims ge 3 ) then begin
   NCDF_DIMINQ, fileID, 2, dimstr, dimsize
  ;help, dimstr
   if (dimstr eq 'time') then begin
     NT=dimsize
   endif else begin
     NZ=dimsize
   endelse
  endif else begin
    if ( finq.ndims ge 4 ) then begin
     NCDF_DIMINQ, fileID, 3, dimstr, dimsize
    ;help, dimstr
     NT=dimsize
    endif
  endelse
 ;help, NX, NY, NZ, NT

  if (keyword_set(LEVEL)) then begin
   if (NZ ge 1) then begin
    varID=ncdf_varid(fileID,'lev')
    ncdf_varget,fileID,varID,levs
   endif
  endif

  if (myVar eq 'WSPD50M') then begin

     varID=ncdf_varid(fileID,'U50M')
     ncdf_varget,fileID,varID,u
     varID=ncdf_varid(fileID,'V50M')
     ncdf_varget,fileID,varID,v
     cdataM = SQRT(u*U + v*v)

  endif else begin

  myVars=STRSPLIT(myVar,'*', /extract)
  s=SIZE(myVars)
  if (s(1) eq 1) then begin
     varID=ncdf_varid(fileID,myVar)
     if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataM,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),TIME]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
     endif else begin
     ;help, TIME
      ncdf_varget,fileID,varID,cdataM,COUNT=[NX,NY,1], OFFSET=[0,0,TIME]
     endelse
  endif else begin
     varID=ncdf_varid(fileID,myVars(0))
     if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataA,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),0]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
     endif else begin
      ncdf_varget,fileID,varID,cdataA
     endelse
     isUNDEF = WHERE(cdataA eq undef)
     if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataA)) then stop
     if (isUNDEF(0) ne -1) then cdataA(isUNDEF) = 0.0 ;!VALUES.F_NAN
     varID=ncdf_varid(fileID,myVars(1))
     if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataB,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),0]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
     endif else begin
      ncdf_varget,fileID,varID,cdataB
     endelse
     isUNDEF = WHERE(cdataB eq undef)
     if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataB)) then stop
     if (isUNDEF(0) ne -1) then cdataB(isUNDEF) = 0.0 ;!VALUES.F_NAN
     cdataM=cdataA*cdataB
     cdataA=0
     cdataB=0
  endelse

  endelse

  if (keyword_set(TRMM_DATA)) then begin
    nLons = 1440
    nLats = 720
    data = FLTARR(nLons,nLats)
    data(*,*) = !VALUES.F_NAN
    data(*,160:559) = cdataM
    cdataM=data
  endif

  if (keyword_set(CMORPH_DATA)) then begin
    nLons = 1440
    nLats = 720
    data = FLTARR(nLons,nLats)
    data(*,*) = !VALUES.F_NAN
    data(*,120:599) = REVERSE(cdataM,2)
    cdataM=data
    isUNDEF = WHERE(cdataM eq MIN(cdataM,/NAN))
    if (isUNDEF(0) ne -1) then cdataM(isUNDEF) = !VALUES.F_NAN
    flip_hemispheres, cdataM
  endif

  isUNDEF = WHERE(cdataM eq undef)
  if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataM)) then stop
  if (isUNDEF(0) ne -1) then cdataM(isUNDEF) = !VALUES.F_NAN

  if (keyword_set(SCALE_FACTOR)) then begin
      print, 'SCALING by: ', SCALE_FACTOR
      cdataM = cdataM*SCALE_FACTOR
  endif

; Check if data is on cubed-sphere grid
  isCubedSphere=0
  if (NY eq 6*NX) then isCubedSphere=1
  if (ARRAY_EQUAL(cdataM,0.0)) then stop
 ;isNAN = WHERE(FINITE(cdataM))
 ;if (N_ELEMENTS(isNAN) ne N_ELEMENTS(cdataM)) then stop 
 ;isUNDEF = WHERE(cdataM eq undef)
 ;if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataM)) then stop
  if (~keyword_set(INFILEID)) then ncdf_close,fileID

  if (isCubedSphere AND keyword_set(regridCubeToLatlon)) then begin
  data=FLTARR(gLons,gLats)
  get_domain, data, LonBeg, LonEnd, LatBeg, LatEnd
  m=size(data)
  cLons=m(1)
  cLats=m(2)
  endif else begin
  cLons=gLons
  cLats=gLats
  endelse

 ;help, isCubedSphere, regridCubeToLatlon
  if (isCubedSphere AND keyword_set(regridCubeToLatlon)) then begin
     print, myVar
     print, min(cdataM,/NAN), max(cdataM,/NAN)
     if (keyword_set(regridResolution)) then begin
        print, 'Reducing Cubed-Sphere Resolution to: ', regridResolution
        print, NX, NY
        t1 = FLTARR(NX,NX)
        t2 = FLTARR(NX,NX)
        t3 = FLTARR(NX,NX)
        t4 = FLTARR(NX,NX)
        t5 = FLTARR(NX,NX)
        t6 = FLTARR(NX,NX)
        for n=0,5 do begin
          i1=0
          i2=NX-1
          j1=NX*n
          j2=j1+NX-1
          if (n eq 0) then t1(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 1) then t2(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 2) then t3(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 3) then t4(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 4) then t5(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 5) then t6(*,*) = cdataM(i1:i2,j1:j2)
        endfor
        NX2 = regridResolution
        XLocs = NX*FINDGEN(NX2)/FLOAT(NX2-1)
        t1 = INTERPOLATE(t1, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t2 = INTERPOLATE(t2, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t3 = INTERPOLATE(t3, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t4 = INTERPOLATE(t4, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t5 = INTERPOLATE(t5, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t6 = INTERPOLATE(t6, XLocs, XLocs, /GRID);, CUBIC=-0.5)
      ; t1 = BILINEAR(t1, XLocs, XLocs)
      ; t2 = BILINEAR(t2, XLocs, XLocs)
      ; t3 = BILINEAR(t3, XLocs, XLocs)
      ; t4 = BILINEAR(t4, XLocs, XLocs)
      ; t5 = BILINEAR(t5, XLocs, XLocs)
      ; t6 = BILINEAR(t6, XLocs, XLocs)
        cdataM=[[t1],[t2],[t3],[t4],[t5],[t6]]
       ;help, cdataM
        NX = NX2
        NY = NX2*6
        print, NX, NY
     endif
if (keyword_set(CONSERV)) then begin
     print, 'Conservative Regridding Cube to Latlon'
     NLONS=NX*4
     NLATS=NX*2 + 1
     sNX=STRING(NX,FORMAT='(i4.4)')
     sNLONS=STRING(NLONS,FORMAT='(i4.4)')
     sNLATS=STRING(NLATS,FORMAT='(i4.4)')
     TILE_FILE='/discover/nobackup/ltakacs/bcs/Ganymed-4_0/Ganymed-4_0_Ostia/Shared/DC'+sNLONS+'xPC'+sNLATS+'_CF'+sNX+'x6C.bin'
     print, TILE_FILE
     spawn, 'ls -l '+TILE_FILE
     openr, tile_lun, TILE_FILE, /GET_LUN, /F77_UNFORMATTED
     NT=0L
     NGRIDS=0
     GRIDNAME=''
     IM=0L
     JM=0L
     readu, tile_lun, NT
     readu, tile_lun, NGRIDS
     help, NT, NGRIDS
     for n=1,NGRIDS do begin
       readu, tile_lun, GRIDNAME
       readu, tile_lun, IM
       readu, tile_lun, JM
       help, IM, JM
     endfor
     IIo=LONARR(NT)
     JJo=LONARR(NT)
      Wo=FLTARR(NT)
     IIi=LONARR(NT)
     JJi=LONARR(NT)
      Wi=FLTARR(NT)
     dummy=''
     readu, tile_lun, dummy ;TYPE
     readu, tile_lun, dummy ;X
     readu, tile_lun, dummy ;Y
   ; Output Grid
     dLON   = 360.0/FLOAT(NLONS)
     dLAT   = 180.0/FLOAT(NLATS-1)
     lonsO  = (-180.0 + 0.5*dLON + FINDGEN(NLONS)*dLON)
     latsO  = ( -90.0 + FINDGEN(NLATS)*dLAT)
     readu, tile_lun,  Wo
                   IIo=Wo
     readu, tile_lun,  Wo
                   JJo=Wo
     readu, tile_lun,  Wo
   ; Input Grid
     readu, tile_lun,  Wi
                   IIi=Wi
     readu, tile_lun,  Wi
                   JJi=Wi
     readu, tile_lun,  Wi
   ; DONE
     close, tile_lun
     free_lun, tile_lun
     help, IIo, JJo, Wo
     print, MIN(IIo), MAX(IIo)
     print, MIN(JJo), MAX(JJo)
     print, MIN( Wo), MAX( Wo)
     help, IIi, JJi, Wi
     print, MIN(IIi), MAX(Iii)
     print, MIN(JJi), MAX(JJi)
     print, MIN( Wi), MAX( Wi)
     cdataO=FLTARR(NLONS,NLATS)
     cdataO(*,*) = 0.0
         FF=FLTARR(NLONS,NLATS)
         FF(*,*) = 0.0
     for N=0,NT-1 do begin
        IO=IIo(N)-1
        JO=JJo(N)-1
       ; Check if I'm in plot Range
        if ( (lonsO(IO) ge LonBeg) and (lonsO(IO) le LonEnd) AND $
             (latsO(JO) ge LatBeg) and (latsO(JO) le LatEnd) ) then begin
         II=IIi(N)-1
         JI=JJi(N)-1
         VALI = cdataM(II,JI)
         if (VALI ne undef) then begin
           W = Wo(N)
           cdataO(IO,JO) = cdataO(IO,JO) + W * VALI
               FF(IO,JO) =     FF(IO,JO) + W
         endif
        endif
     endfor
     izeros = WHERE(FF eq 0.0, COMPLEMENT=isVALID)
     cdataO(isVALID)=cdataO(isVALID)/FF(isVALID)
     cdataO(izeros) =!VALUES.F_NAN
     cdataM=cdataO
     isCubedSphere=0
     PIXEL=1
     cLons=gLons
     cLats=gLats
endif else begin
     print, 'Regridding Cube to Latlon'
;  Fix missing data
     imiss = WHERE(cdataM eq undef)
     if (imiss(0) ne -1) then cdataM(imiss) = !VALUES.F_NAN

; Setup the mapping
     if (~keyword_set(PROJECT)) then begin
      SET_PLOT, 'Z'
      DEVICE, SET_RESOLUTION=[cLons, cLats]
      loadct, 0
      erase, 0
      maplimit=STRING(LatBeg)+STRING(LonBeg)+STRING(LatEnd)+STRING(LonEnd)
      strMap = "MAP_SET, XMARGIN=0, YMARGIN=0, /NOBORDER, LIMIT=["+STRJOIN(STRSPLIT(maplimit,/EXTRACT),',')+"]"
      print, strMap
      res=EXECUTE(strMap)
     endif
     bdataM=BYTSCL(cdataM)
     if (keyword_set(DATA_LEVELS)) then begin
       dlevs=DATA_LEVELS
       nlevs=N_ELEMENTS(dlevs)
       blevs=254L*(FINDGEN(nlevs)/(nlevs-1))+1
      ;help, dlevs, blevs
     endif else begin
       if (keyword_set(RANGE)) then begin
         cMIN=RANGE[0]
         cMAX=RANGE[1]
       endif else begin
         cMIN=MIN(cdataM, /NAN)
         cMAX=MAX(cdataM, /NAN)
       endelse
      ;help, cMIN, cMAX
       dlevs=(FINDGEN(255)/254.0)*(cMAX-cMIN) + cMIN
       if (keyword_set(LOG)) then begin
          logBins=cMIN + [0.345, 0.402, 0.468, 0.545, 0.634, 0.7381, 0.859129, 1, 1.16397, 1.35483, 1.57698, 1.83556, 2.13654, 2.48687, 2.89464, 3.36928, 3.92175, 4.5648, 5.31329, 6.18452, 7.1986, 8.37896, 9.75287, 11.3521, 13.2135, 15.3801, 17.902, 20.8374, 24.2541, 28.2311, 32.8602, 38.2483, 44.5199, 51.8198, 60.3168, 70.207, 81.7189, 95.1185, 110.715, 128.869, 150.0, 174.6, 203.2]*(cMAX-cMIN)/203.2
          XLocs = (N_ELEMENTS(logBins)-1)*FINDGEN(255)/FLOAT(255-1)
          dlevs = INTERPOLATE(logBins, XLocs, /GRID)
          print, dlevs
       endif
       blevs=(FINDGEN(255)+1)
     endelse
     NCOLORS=255
     if (keyword_set(NCOLORS)) then begin
        XLocs = (N_ELEMENTS(dlevs)-1)*FINDGEN(NCOLORS)/FLOAT(NCOLORS-1)
        dlevs = INTERPOLATE(dlevs, XLocs)
        nIncr = FLOAT(N_ELEMENTS(dlevs)-1)/3.0
        dlevs = SMOOTH(dlevs,nIncr)
        blevs=(FINDGEN(NCOLORS)+1)
     endif
     print, dlevs

; map the data
if (1) then begin
     iplot=0
   ; CUBE FACES
     for tile=1,6 do begin

     inDomain=0
     print, 'Checking for data in Cube Tile '+STRTRIM(STRING(tile),1)
     gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20090716/CUBED_GRIDS/c48_grid.tile'+STRTRIM(STRING(tile),1)+'.nc'
     if (keyword_set(STRETCH)) then gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20160301/CUBED_GRIDS/c48_grid_stretch.tile'+STRTRIM(STRING(tile),1)+'.nc'
     fileID = ncdf_open(gridfile, /NOWRITE)
     varID=ncdf_varid(fileID,'x')
     ncdf_varget,fileID,varID,Xverts
     iNeg=WHERE(Xverts gt 180)
     if (iNeg(0) ne -1) then Xverts(iNeg)=Xverts(iNeg)-360
     varID=ncdf_varid(fileID,'y')
     ncdf_varget,fileID,varID,Yverts
     ncdf_close, fileID
     ifind = WHERE( (Xverts ge LonBeg) AND (Xverts le LonEnd) )
     if (ifind(0) ne -1) then inDomain=1
     ifind = WHERE( (Yverts ge LatBeg) AND (Yverts le LatEnd) )
     if (ifind(0) ne -1) then inDomain=inDomain+1
     print, MIN(Xverts), MAX(Xverts)

     if (inDomain eq 2) then begin
       gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20090716/CUBED_GRIDS/c'+STRTRIM(STRING(NX),1)+'_grid.tile'+STRTRIM(STRING(tile),1)+'.nc'
       if (keyword_set(STRETCH)) then gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20160301/CUBED_GRIDS/c'+STRTRIM(STRING(NX),1)+'_grid_stretch.tile'+STRTRIM(STRING(tile),1)+'.nc'
       fileID = ncdf_open(gridfile, /NOWRITE)
       varID=ncdf_varid(fileID,'x')
       ncdf_varget,fileID,varID,Xverts
       iNeg=WHERE(Xverts gt 180)
       if (iNeg(0) ne -1) then Xverts(iNeg)=Xverts(iNeg)-360
       varID=ncdf_varid(fileID,'y')
       ncdf_varget,fileID,varID,Yverts
     ; help, Xverts, Yverts
       print, MIN(Xverts), MAX(Xverts)
     ; print, MIN(Yverts), MAX(Yverts)
       print, 'Plotting Cube Tiles '+STRTRIM(STRING(tile),1)
         for i = 0, NX-1 do begin
            if ((i mod 288) eq 0) then print, 'Plotting Cube Tiles '+STRTRIM(STRING(tile),1)+' i='+STRTRIM(STRING(i),1)
            for j = 0, NX-1 do begin
              iS=i*2
              jS=j*2
              iC=i
              jC=j + (tile-1)*NX
             ; Check if I'm in plot Range
              if ( (MAX(Xverts(iS:iS+2,jS:jS+2)) ge LonBeg) and (MIN(Xverts(iS:iS+2,jS:jS+2)) le LonEnd) ) then begin
              if ( (MAX(Yverts(iS:iS+2,jS:jS+2)) ge LatBeg) and (MIN(Yverts(iS:iS+2,jS:jS+2)) le LatEnd) ) then begin
               xsym = [Xverts(iS,jS),Xverts(iS,jS+2),Xverts(iS+2,jS+2),Xverts(iS+2,jS)]
               ysym = [Yverts(iS,jS),Yverts(iS,jS+2),Yverts(iS+2,jS+2),Yverts(iS+2,jS)]
               ifind=WHERE(dlevs ge cdataM(iC,jC))
               myColor = ifind(0)+1
               if (FINITE(cdataM(iC,jC)) eq 0) then myColor = 0

               if ( (keyword_set(STRETCH)) AND $
                    ((MIN(Yverts(iS:iS+2,jS:jS+2)) gt 85) OR MAX(Yverts(iS:iS+2,jS:jS+2)) lt -85) ) then begin
            ;    SORTING = SORT(xsym)
            ;    ysym = ysym(SORTING)
            ;    xsym = xsym(SORTING)
                 if (iplot eq 0) then begin
                ;help, iplot
                ;print, xsym, ysym
                ;polyfill, xsym, ysym, /data, color=myColor
                 endif
                 iplot=iplot+1
               endif else begin
                 polyfill, xsym, ysym, /data, color=myColor
               endelse
            ;  plots,  xsym, ysym, /data, color=myColor, thick=3

              endif
              endif
            endfor
         endfor
       ncdf_close, fileID
     endif
     endfor ; tiles
endif else begin
     varID=ncdf_varid(fileID,'CLONS')
     ncdf_varget,fileID,varID,LONS
     varID=ncdf_varid(fileID,'CLATS')
     ncdf_varget,fileID,varID,LATS
     ncdf_close,fileID
     LONS=LONS*(180.0/!DPI)
     LATS=LATS*(180.0/!DPI)
     mySym=SYM(1)
     for i=0,NX-1 do begin
     if ((i mod 288) eq 0) then print, 'Plotting Cube Points i='+STRTRIM(STRING(i),1)
     for j=0,NY-1 do begin
        ifind=WHERE(dlevs ge cdataM(i,j)) + 1
       ;plots,LONS[i,j],LATS[i,j],COLOR=ifind(0),PSYM=3
        plots,LONS[i,j],LATS[i,j],COLOR=ifind(0),PSYM=mySym,SYMSIZE=0.5
     endfor
     endfor
endelse
    cdataM=FLOAT(TVRD())
   ;help, cdataM
    if (0) then begin
     for i=0,cLons-1 do begin
     for j=0,cLats-1 do begin
        ifind=WHERE(blevs ge cdataM(i,j))
        cdataM(i,j) = dlevs(ifind(0))
     endfor
     endfor
    endif else begin
      cdataB=cdataM
      for n=1,N_ELEMENTS(dlevs) do begin
        ifind = WHERE(cdataB ge n)
        if (ifind(0) ne -1) then begin
           cdataM(ifind) = dlevs(n-1)
        endif
      endfor
    endelse
    imiss = WHERE(cdataB eq 0)
    if (imiss(0) ne -1) then cdataM(imiss)=!VALUES.F_NAN
    help, imiss
    ifind = WHERE(cdataM lt dlevs(0))
    cdataM(ifind) = dlevs(0)
    print, min(cdataM,/NAN), max(cdataM,/NAN)
  endelse
  endif


endif else begin
  print, myVar
  print, myfile
  stop
endelse
imiss = WHERE(cdataM eq undef)
if (imiss(0) ne -1) then cdataM(imiss) = !VALUES.F_NAN
s= SIZE(cdataM)
is3D=0
if (s(0) eq 3) then is3D=1
is2D=0
if (s(0) eq 2) then is2D=1
is1D=0
if (s(0) eq 1) then is1D=1
nLevs=0
if (is3D) then nLevs=s(3)

if ( (not isCubedSphere) OR (isCubedSphere and keyword_set(regridCubeToLatlon)) ) then begin
if (~keyword_set(PROJECT)) then begin
; Interpolate
s= SIZE(cdataM)
nLons=DOUBLE(s(1))
nLats=DOUBLE(s(2))
if ( (nLons ne cLons) OR (nLats ne cLats) ) then begin
  if (is3D) then begin
     XLocs = nLons*FINDGEN(cLons)/FLOAT(cLons-1)
     YLocs = nLats*FINDGEN(cLats)/FLOAT(cLats-1)
     ZLocs = FINDGEN(nLevs)
     cdataM = INTERPOLATE(cdataM, XLocs, YLocs, ZLocs, /GRID)
  endif else begin
     if (is2D) then begin
       if (keyword_set(PIXEL)) then begin
       print, 'PIXEL Interpolation'
       XLocs = LONG(nLons*FINDGEN(cLons)/FLOAT(cLons-1))
       YLocs = LONG(nLats*FINDGEN(cLats)/FLOAT(cLats-1))
       cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
       endif else begin
       print, 'SMOOTH Interpolation'
      ;if (nLons GT cLons) then begin
      ;  oLons = nLons
      ;  oLats = nLats
      ;  iLons = oLons/2.0
      ;  iLats = (oLats-1)/2.0 + 1.0
      ;  while (iLons ge cLons) do begin  
      ;    help, cdataM
      ;    print, iLons, iLats, cLons, cLats
      ;    XLocs = (oLons*FINDGEN(iLons)/FLOAT(iLons-1))
      ;    YLocs = (oLats*FINDGEN(iLats)/FLOAT(iLats-1))
      ;    cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
      ;    oLons = iLons
      ;    oLats = iLats
      ;    iLons = oLons/2.0
      ;    iLats = (oLats-1)/2.0 + 1.0
      ;  endwhile
      ;endif else begin
         XLocs = (nLons*FINDGEN(cLons)/FLOAT(cLons-1))
         YLocs = (nLats*FINDGEN(cLats)/FLOAT(cLats-1))
         cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
      ;  help, XLocs, YLocs
      ;  help, cdataM
      ;endelse
       endelse
     endif
  endelse
endif
;help, cdataM

; Subdomain
if ( (not isCubedSphere) OR (isCubedSphere and ~keyword_set(regridCubeToLatlon)) ) then begin
if (is3D) then begin
  if (keyword_set(FLIP)) then flip_hemispheres, cdataM
  for k=0,nLevs-1 do begin
     data2D=cdataM(*,*,k)
     get_domain, data2D, LonBeg, LonEnd, LatBeg, LatEnd
     if (k eq 0) then begin
        s=size(data2D)
        sLons=s(1)
        sLats=s(2)
        data3D=FLTARR(sLons,sLats,nLevs)
     endif
    ; Flip the vertical so 0:nLevs-1 goes from surface:top
     data3D(*,*,nLevs-k-1)=data2D
  endfor
  data=cdataM
  data2D=0
  data3D=0
endif else begin
  if (is2D) then begin
     if (keyword_set(FLIP)) then flip_hemispheres, cdataM
     get_domain, cdataM, LonBeg, LonEnd, LatBeg, LatEnd
     data=cdataM
  endif else begin
   if (myVar eq "lon") then begin
      dlon = FLOAT(LonEnd-LonBeg)/FLOAT(cLons)
      data = LonBeg + dlon*FINDGEN(cLons)
   endif
   if (myVar eq "lat") then begin
      dlat = FLOAT(LatEnd-LatBeg)/FLOAT(cLats)
      data = LatBeg + dlat*FINDGEN(cLats)
   endif
  endelse
endelse
cdataM=0
endif else begin
  data=cDataM
endelse

endif else begin
  data=cDataM
endelse

endif else begin ; PROJECT
  data=cDataM
endelse

print, myVar
help, data
print, min(data, /NAN), max(data, /NAN)

end

; cube_to_latlon.pro
pro cube_to_latlon

set_plot, 'Z'

LonBeg = -180
LonEnd =  180
LatBeg =  -90
LatEnd =   90
lonCen = 0.5*(LonBeg+LonEnd)
LatCen = 0.0
undef = 1.e15
pngImgIdim=1440*4.0
pngImgJdim=pngImgIdim*0.5 + 1
boundaryImageObj = get_map_image_new(LatBeg, LatEnd, LonBeg, LonEnd, LonCen, LatCen, $
                   /NOBORDER, HORIZ_DIMS=[pngImgIdim, pngImgJdim])
elevfile='/discover/nobackup/projects/gmao/osse2/stage/c1440_NR/DATA/0.0625_deg/const/const_2d_asm_Nx/Y2005/M11/D01/c1440_NR.const_2d_asm_Nx.20051101.nc4'
read_and_interpolate, elevfile, 'lon', lons, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef
read_and_interpolate, elevfile, 'lat', lats, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef
myfile = "/discover/nobackup/projects/gmao/osse2/G5FCST/forecast-20121024_150000-C3072-4390102/G5FCST.inst1_2d_asm_Nx.20121026_1200z.nc4"
;myfile = "/discover/nobackup/projects/gmao/osse2/c5760_NR/scratch-3550451.out/c5760_NR.inst10mn_2d_met1_Nx.20120616_0840z.nc4"
;myfile = "/discover/nobackup/projects/gmao/osse2/H10UC90R72/scratch-3861154/H10UC90R72.inst30mn_2d_met1_Nx.20150215_0000z.nc4"
read_and_interpolate_cube, myfile, 'U250', u250, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef, /regridCubeToLatlon
read_and_interpolate_cube, myfile, 'V250', v250, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef, /regridCubeToLatlon

; write netcdf file
; Create a new NetCDF file with the filename inquire.nc:
id = NCDF_CREATE('c3072_NR_U250-V250.nc', /CLOBBER)
; Fill the file with default values:
NCDF_CONTROL, id, /FILL
xid = NCDF_DIMDEF(id, 'x', pngImgIdim)
yid = NCDF_DIMDEF(id, 'y', pngImgJdim)
tid = NCDF_DIMDEF(id, 't', 1)
; Define variables:
latid  = NCDF_VARDEF(id,  'lat', [yid], /FLOAT)
lonid  = NCDF_VARDEF(id,  'lon', [xid], /FLOAT)
timeid = NCDF_VARDEF(id, 'time', [tid], /FLOAT)
uid    = NCDF_VARDEF(id, 'U250', [xid,yid], /FLOAT)
vid    = NCDF_VARDEF(id, 'V250', [xid,yid], /FLOAT)
NCDF_CONTROL, id, /ENDEF
; Input data:
NCDF_VARPUT, id, timeid, 0
NCDF_VARPUT, id,  lonid, lons
NCDF_VARPUT, id,  latid, lats
NCDF_VARPUT, id,    uid, u250
NCDF_VARPUT, id,    vid, v250
NCDF_CLOSE, id ; Close the NetCDF file.

end

GrADS

To use GrADS to display cube-sphere grid data in a NetCDF file, the path should be in the PATH environment: /discover/nobackup/projects/gmao/share/dasilva/opengrads/Contents/grads

Then the users can follow the steps. First, we assume a version of GEOS-5 is available.

1. run $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/configure. This will produce a .quickplotrc in that directory
2. source $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/.quickplotrc
3. run grads
4. sdfopen file.nc4
5. run command dc to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument.