PRO Geostrophic ; =============================================================================== ; Program to plot pressure and wind vector composite fields ; =============================================================================== ; ---------------------------------------------------------------------- ; Set-up the colours and the postscript output ; ---------------------------------------------------------------------- PRINT, 'Setting-up the colours and the postscript device' r = [0,255,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184] g = [0,255,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0] b = [0,255,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186] TVLCT, r,g,b !P.POSITION = [0.15, 0.15, 0.95, 0.95] SET_PLOT, 'ps' PostscriptFile = 'Geostrophic.ps' DEVICE, FILE=PostscriptFile, /PORTRAIT, /COLOR, XSIZE=40.0, YSIZE=25.0, $ XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8, /ENCAPSULATED plotheight = 5500.0 ; This is the height to look at ; ---------------------------------------------------------------------- ; Read-in the wind, pressure and orography data ; ---------------------------------------------------------------------- PRINT, 'Reading in the data' fileid = NCDF_OPEN('Fc.nc') varidx = NCDF_VARID(fileid, 'longitude_1') varidy = NCDF_VARID(fileid, 'latitude') varidz = NCDF_VARID(fileid, 'hybrid_ht') varidu = NCDF_VARID(fileid, 'u') varidv = NCDF_VARID(fileid, 'v') varidp = NCDF_VARID(fileid, 'field7') ; pressure varidorog = NCDF_VARID(fileid, 'ht') ; orography NCDF_VARGET, fileid, varidx, x_deg NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varidz, vertlevs NCDF_VARGET, fileid, varidu, uwind_all NCDF_VARGET, fileid, varidv, vwind_all NCDF_VARGET, fileid, varidp, press_all NCDF_VARGET, fileid, varidorog, orog NCDF_CLOSE, fileid ; ---------------------------------------------------------------------- ; Calculate the height field (fn of long and lat) from sea level ; ---------------------------------------------------------------------- PRINT, 'Calculating the height field' eta = FLTARR(38) ztop = vertlevs(37) eta(0:37) = vertlevs(0:37) / ztop interface = 30-1 height = FLTARR(288+1,217,38) FOR x = 0, 288-1 DO BEGIN FOR y = 0, 217-1 DO BEGIN height(x,y,0:37) = eta(0:37) * ztop height(x,y,0:interface) = height(x,y,0:interface) + $ (1.0 - eta(0:interface) / eta(interface)) * $ (1.0 - eta(0:interface) / eta(interface)) * $ orog(x,y) ENDFOR ENDFOR ; ---------------------------------------------------------------------- ; Set-up arrays and interpolate the data to the specified height ; ---------------------------------------------------------------------- PRINT, 'Interpolating the data' uwind = FLTARR(288+1,217) vwind = FLTARR(288+1,217) press = FLTARR(288+1,217) result = 0.0 FOR x = 0, 288-1 DO BEGIN FOR y = 0, 216-1 DO BEGIN Find, height, uwind_all, x, y, plotheight, result uwind(x,y) = result Find, height, vwind_all, x, y, plotheight, result vwind(x,y) = result Find, height, press_all, x, y, plotheight, result press(x,y) = result ENDFOR ENDFOR ; ---------------------------------------------------------------------- ; Wrap the data around the Greenwich meridian ; ---------------------------------------------------------------------- PRINT, 'Wrapping the data over the Greenwich meridian' x_degrees = FLTARR(288+1) x_degrees(0:288-1) = x_deg(0:288-1) x_degrees(288) = x_degrees(287) + x_degrees(287) - x_degrees(286) uwind(288,*) = uwind(0,*) vwind(288,*) = vwind(0,*) press(288,*) = press(0,*) ; ---------------------------------------------------------------------- ; Set-up the map projection ; ---------------------------------------------------------------------- PRINT, 'Setting-up the map projection' minx = -60.0 maxx = 30.0 miny = 20.0 maxy = 70.0 MAP_SET, 0.0, 0.0, LIMIT=[maxy,minx,miny,maxx], /CONTINENTS, /CYLINDRICAL, /NOBORDER ; Midlats ; ---------------------------------------------------------------------- ; Plot the pressure data - colours ; ---------------------------------------------------------------------- ;nlevs = 16 ;min = MIN(press(*,*)) ;max = MAX(press(*,*)) ;levs = min + FINDGEN(nlevs) * (max-min) / nlevs ;CONTOUR, press(*,*), x_degrees, y_degrees, LEVELS=levs, XSTYLE=5, YSTYLE=5, /CELL_FILL, C_COLORS=INDGEN(15)+2, /OVERPLOT ; ---------------------------------------------------------------------- ; Plot the pressure data - contours ; ---------------------------------------------------------------------- PRINT, 'Plotting the pressure data' nlevs = 24 min = MIN(press(0:288,0:215)) max = MAX(press(0:288,0:215)) levs = min + FINDGEN(nlevs) * (max-min) / nlevs PRINT, ' Min pressure = ', min PRINT, ' Max pressure = ', max labs = STRTRIM(STRING(levs),2) CONTOUR, press(*,*), x_degrees, y_degrees, LEVELS=levs, XSTYLE=5, YSTYLE=5, C_ANNOTATION=labs, COLOR=13, /OVERPLOT ; ---------------------------------------------------------------------- ; Plot the wind arrows ; ---------------------------------------------------------------------- PRINT, 'Plotting the wind data' VELOVECT, uwind(0:288,0:215), vwind(0:288,0:215), x_degrees(0:288), y_degrees(0:215), LENGTH=2, THICK=1.0, /OVERPLOT, /NOERASE ;----------------------------------------------------------- ; Add continents ;----------------------------------------------------------- PRINT, 'Adding the continents' MAP_CONTINENTS, /NOERASE ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- PRINT, 'Adding the axes' blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE='degrees' AXIS, COLOR=0, XAXIS=1, TICKLEN=-0.015, XRANGE=[minx,maxx], XTICKNAME=blank AXIS, COLOR=0, YAXIS=0, TICKLEN=-0.015, YRANGE=[miny,maxy], YTITLE='degrees' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank DEVICE, /CLOSE END PRO Find, HeightArr, ValueArr, xx, yy, height, value lev = 0 found = 0 WHILE (lev LT 37) && (found EQ 0) DO BEGIN ; PRINT, lev, height, HeightArr(xx, yy, lev), HeightArr(xx, yy, lev+1) IF (height GT HeightArr(xx, yy, lev)) && (height LE HeightArr(xx, yy, lev+1)) THEN BEGIN found = 1 ENDIF ELSE BEGIN lev = lev + 1 ENDELSE ENDWHILE IF (found EQ 1) THEN BEGIN ; Do linear interpolation alpha = (height - HeightArr(xx, yy, lev)) / (HeightArr(xx, yy, lev+1) - HeightArr(xx, yy, lev)) beta = (HeightArr(xx, yy, lev+1) - height) / (HeightArr(xx, yy, lev+1) - HeightArr(xx, yy, lev)) value = alpha * ValueArr(xx, yy, lev+1) + beta * ValueArr(xx, yy, lev) ENDIF ELSE BEGIN value = ValueArr(xx, yy, 0) ;; Use surface values PRINT, 'Warning - ', xx, yy ENDELSE END