;Developed by Wayne Feltz 3/15/95 ;Space Science and Engineering Center (SSEC) University of Wisconsin ; This idl program was developed to read aeri binary retrieval format (brf) ;files and produce atmospheric time cross sections. These cross sections may ;be temperature, potential temperature, virtual temperature, mixing ratio, and ;relative humidity. Total precipitable water is also calculated and stored in ;an array. This can be ploted at the idl prompt after aerixsec has been ;terminated. (Useful for comparing to microwave radiometer data) ;If any questions please contact Wayne Feltz, waynef@ssec.wisc.edu or ;608-265-6283 ;The following modules are used with this program: ;filter,hourform,indaxe,relhum,satmix,svpice,svpwat,totpre,virtemp,vline2,wvmr ;*************************************************************************** ;Works with no ozone brf file 7/7/95 WFF ;*************************************************************************** !P.font=-1 COMMON Share1, p, time close, 1 set_plot, 'X' !P.multi=[0,1,2,0,0] !x.margin = [8,6] !y.margin = [4,4] window, 0, colors=256, xsize = 700, ysize = 700,retain=2 ;average span cin=255 ;read data as1=5 as2=5 filename=pickfile(path = '/home3/waynef/camexiii/andros/cdffile/', filter='98*pr.cdf') ;############################################################# ; READ IN NETCDF AERI PHYSICAL RETRIEVALS fid = ncdf_open(filename) ncdf_varget, fid, 'height', alt ncdf_varget, fid, 'pressure', pr ncdf_varget, fid, 'ambientTemp', temp ncdf_varget, fid, 'dewpointTemp', dtemp ncdf_varget, fid, 'AERItime', time ncdf_varget, fid, 'cloudBaseHeight', cbh ncdf_varget, fid, 'Date', date stime=string(time) recno=n_elements(time) sdate=strcompress(string(date),/remove_all) recst=3 pr=reverse(pr(0:19,*)) temp=reverse(temp(0:19,*)) dtemp=reverse(dtemp(0:19,*)) mixr=fltarr(20,recno) alt=reverse(alt) ;############################################################# for i = 0, recno - 1 do mixr(*,i)=wvmr(pr(*,i), $ temp(*,i),dtemp(*,i)) minalt = 0.0 maxalt = max(alt(where(temp(0:19,0) GT 0.0))) top=2250. bot=0.0 ibot=19 itop=where(alt eq top) itop = fix(fix(itop),0) height = 20 - where(alt eq top) height = fix(fix(height),0) ;matrix to contour surfta = fltarr(recno, height) surftha = fltarr(recno, height) surfthae = fltarr(recno, height) surfma = fltarr(recno, height) surftva = fltarr(recno, height) surfrha = fltarr(recno, height) surfva_ns= fltarr(recno, height) ;first record p = pr(itop:ibot,recst) al=alt(itop:ibot) surfta(recst, *) = temp(itop:ibot,recst) surfma(recst,*) = mixr(itop:ibot,recst) ;Interpolate AERI temperature and mixing ratio data for i = recst+1, recno-1 do begin surfta(i, *) = temp(itop:ibot,i) surfma(i,*) = mixr(itop:ibot,i) endfor result=where(surfta eq -8888.00,count) if (count gt 1) then surfta(result)=-500.0 if (count gt 1) then surfma(result)=-500.0 for i = recst, recno-1 do surfva_ns(i,*)=virtemp(surfta(i,*),surfma(i,*))-273.16 ;smooth surf with a running average filter for i = 0, height-1 do $ surfta(recst:recno-1,i) = filter(surfta(recst:recno-1, i), as1, recno-recst) for i = 0, height-1 do $ surfma(recst:recno-1,i) = filter(surfma(recst:recno-1, i), as2, recno-recst) ;Calculate total preciptable water for later use ;Calculate theta for AERI cross section preh2o=fltarr(recno) pbot = (max(pr(*,0)))-1 for i = recst, recno-1 do begin totpre,p, surfma(i,*), height, aeri_ini_temp,aeri_total_temp surftha(i,*)=surfta(i,*)*(pbot/p)^.286 surftva(i,*)=virtemp(surfta(i,*),surfma(i,*)) surfrha(i,*)=relhum(p,surfta(i,*),surfma(i,*)) ; surfthae(i,*)=thetae(p,surfta(i,*),surfma(i,*),pbot) preh2o(i)=aeri_total_temp endfor result2=where((surfrha lt 0.0) or (surfrha gt 100.0),count) if (count gt 1) then surfrha(result2)=-500.0 ;pressure axes (2) altl = [0.0, 250.0, 500.0, 750.0, 1000.0, 1250.0, 1500.0, 1750.0, 2000.0, $ 2250.0] nintery = (size(altl))(1) ytname = make_array(nintery, /string) FOR i = 0, nintery-1 DO $ ytname(i) = string(fix(altl(i)), format = "(i)") ;desired time ticks (2) st = 2*(fix(min(time(recst:recno-1))/2.0)+1) et = 2*(fix(max(time(recst:recno-1))/2.0)) noet = ((et-st)/2.0)+1 noat = noet + 2 h = fltarr(noat) h(0)=min(time(recst:recno-1)) h(1)=st for i= 2, noat-2 do h(i)=h(i-1)+2.0 h(noat-1)=max(time(recst:recno-1)) ninterx = (size(h))(1) ;time tick name xtname = make_array(ninterx, /string) FOR i = 0, ninterx-1 DO $ xtname(i)=string(fix(hourform(h(i))), format="(i)") ;index of the desired time tick in the raw x range [0,nprof-1] indx = indaxe(time, h, ninterx, 1) winnum=0 repeat begin winnum=winnum+1 ; Pick variable field to be contoured field=1 print, 'Temperature (1)' print, 'Potential Temperature (2)' print, 'Mixing Ratio (3)' print, 'Virtual Temperature (4)' print, 'Relative Humidity (5)' repeat begin read, prompt='Please pick first field to be contoured:',field1 read, prompt='Please pick second field to be contoured:',field2 endrep until (field ge 1) and (field le 5) case field1 of 1:Begin xsec1=surfta type1='!17Wisconsin AERI Sounding System !C!DTEMPERATURE (K) ' label='t' End 2:Begin xsec1=surftha type1='!17Wisconsin AERI Sounding System !C!DPOTENTIAL TEMPERATURE (K) ' label='t' End 3:Begin xsec1=surfma type1='!17Wisconsin AERI Sounding System !C!DMIXING RATIO (G/KG) ' minl1=0 dc1=1.0 ncont1=24 label='m' End 4:Begin xsec1=surftva type1='!17Wisconsin AERI Sounding System !C!DVIRTUAL TEMPERATURE (K) ' label='t' End 5:Begin xsec1=surfrha type1='!17Wisconsin AERI Sounding System !C!DRelative Humidity (%) ' minl1=10 dc1=5 ncont1=18 label='m' End endcase case field2 of 1:Begin xsec2=surfta type2='!17Wisconsin AERI Sounding System !C!DTEMPERATURE (K) ' End 2:Begin xsec2=surftha type2='!17Wisconsin AERI Sounding System !C!DPOTENTIAL TEMPERATURE (K) ' End 3:Begin xsec2=surfma type2='!17Wisconsin AERI Sounding System !C!DMIXING RATIO (G/KG) ' minl2=0.0 dc2=1.0 ncont2=24 End 4:Begin xsec2=surftva type2='!17Wisconsin AERI Sounding System !C!DVIRTUAL TEMPERATURE (K) ' End 5:Begin xsec2=surfrha type2='!17Wisconsin AERI Sounding System !C!DRelative Humidity (%) ' minl2=5.0 dc2=5 ncont2=19 End endcase ;contour levels if (field1 eq 1) or (field1 eq 2) or (field1 eq 4) then begin minl1 = 280 dc1 = 1.0 ncont1=30 endif col = reverse(indgen(ncont1)*255/ncont1) levels1 = findgen(ncont1)*dc1 + float(minl1) ;contour levels ;ncont = 5 if (field2 eq 1) or (field2 eq 2) or (field2 eq 4) then begin ; minl2 = min(xsec2(recst:recno-1, zmax:zmin)) minl2 = 290 ncont2=30 dc2 = 1.0 endif col = reverse(indgen(ncont2)*255/ncont2) levels2 = findgen(ncont2)*dc2 + float(minl2) ;load color table loadct, 33 tvlct,r,b,g,/get r(cin)=255 g(cin)=255 b(cin)=255 tvlct,r,b,g ; blank out data between missing times for j=0,recno-2 do begin dift = time(j+1)-time(j) if dift gt .75 then begin xsec1(j:j+1,*)=-500 xsec2(j:j+1,*)=-500 endif endfor read, prompt = 'Display on the screen (0) or into the Postscript file "idl.ps" (1): ', kind ;display the contour IF kind NE 0 THEN BEGIN read, prompt = 'Landscape (1) or Portrait (2): ', kind IF kind EQ 1 THEN BEGIN set_plot,'ps' device,/landscape,/color ENDIF ELSE BEGIN set_plot,'ps' device,/portrait,/color endelse ENDIF contour, xsec1(recst:recno-1, *), levels = levels1, /follow, /fill,xticks=ninterx-1, $ xtickv = h, xrange = [time(recst), time(recno-1)], $ xstyle = 1, yticks = nintery-1, xtickname=xtname, $ ystyle = 1, ytickname = ytname, ytickv = altl, $ yrange = [altl(0), altl(nintery-1)], charsize=1.2,$ background = cin, color = 0, xtitle = 'Time (UTC)', $ title = type1 + ' AERI CAMEX III ANDROS ISLAND '+ sdate , $ xgridstyle = 1, xticklen = 1.0, ygridstyle = 1, yticklen = 1.0, thick = 1.0, $ ytitle = 'Altitude (m)', font = 1, time(recst:recno-1), al contour, xsec1(recst:recno-1, *), levels = levels1, /follow, /overplot, $ time(recst:recno-1), al, thick = 2.0, c_colors=0 contour, xsec2(recst:recno-1, *), levels = levels2, /follow, /fill,xticks=ninterx-1, $ xtickv = h, xrange = [time(recst), time(recno-1)], $ xstyle = 1, yticks = nintery-1, xtickname=xtname,$ ystyle = 1, ytickname = ytname, ytickv = altl, $ yrange = [altl(0), altl(nintery-1)], charsize=1.2,$ background = cin, color = 0, xtitle = 'Time (UTC)', $ title = type2 + ' AERI CAMEX III ANDROS ISLAND '+ sdate, $ xgridstyle = 1, xticklen = 1.0, ygridstyle = 1, yticklen = 1.0, thick = 1.0, $ ytitle = 'Altitude (m)', font = 1, time(recst:recno-1), al contour, xsec2(recst:recno-1, *), levels = levels2, /follow, /overplot, $ time(recst:recno-1), al, thick = 2.0, c_colors=0 if (kind lt 1) then begin ;repeat begin ; print, 'Would you like to output the window into a graphics file ' ; read, prompt='yes=0 no=1?', gif ;endrep until (gif eq 0) or (gif eq 1) gif = 0 ;Output GIF or BMP files if needed if gif eq 0 then begin print, 'GIF = 0' print, 'BMP = 1' type=0 ; repeat begin ; read, prompt='Please decide on type of file output:', type ; endrep until (type ge 0) and (type le 1) case type of 0:begin tvlct,r,g,b,/Get giffile=tvrd() write_gif, '~/public_html/camexiii/archive/quicklooks/'+sdate+'/'+sdate+label+'.gif', giffile, r, g, b end 1:begin tvlct,r,g,b,/Get bmpfile=tvrd() write_bmp, 'idl.bmp', bmpfile, r, g, b end endcase endif endif IF kind GT 0 THEN device,/close repeat begin read, prompt='Would you like to display two new fields (0=yes) (1=no):',disp endrep until (disp ge 0) and (disp le 1) endrep until (disp eq 1) ncdf_close,fid end