; 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