;--------------------------------------------------------------------- ; ; This program generates a skew-T, log p thermodynamic diagram. This ; program was derived to reproduce the USAF skew-T, log p diagram ; (form DOD-WPC 9-16-1 current as of March 1978). ; ;--------------------------------------------------------------------- pro skewt_grid_fg,position_sonde,pydim,p0,xmin,xmax,ymin,ymax ; These are the functions needed by this subroutine forward_function skewty,skewtx,tda,satlft trange=[-40., 40.] ;temperature range prange=[1050., 100] ;pressure range title='Skew-T, Log(P) Diagram' ; --- Define absoulute x,y max/min bounds corresponding to the outer ; --- edges of the diagram. These are computed by inverting the ; --- appropriate pressures and temperatures at the corners of the diagram. ymax=skewty(prange[0]) ;bottom (1050) ymin=skewty(prange[1]) ;top (100) xmin=skewtx(trange[0],skewty(prange[0])) ;left (-40,1050) xmax=skewtx(trange[1],skewty(prange[0])) ;right (40,1050) ; --- Set up the plot size and position of sounding ; Passing in a position vector and window size if position_sonde ne !NULL then begin x0=position_sonde[0] y0=position_sonde[1] x1=position_sonde[2] y1=position_sonde[3] ; Calculating position vector endif else begin pydim=709 & pxdim=1280 daspect = FLOAT(pydim)/FLOAT(pxdim) * (trange(1)-trange(0))/80.0 margin = 0.1 aspect = 1.0 ; A square x0 = 0.50 - (0.5 - margin)*(daspect/aspect) y0 = margin x1 = 0.50 + (0.75 - margin)*(daspect/aspect) y1 = 1.0 - margin endelse ; Sonde position p_position=[x0, y0, x1, y1] ; Temperature xtitle position temp_title_position=[((x1-x0)/2)+x0,y0-0.07] ;--------------------------------------------------------------------- ; --- DRAW HORIZONTAL ISOBARS., LABEL VERTICAL AXIS ;--------------------------------------------------------------------- ; This draws the ISOBARS p0=plot(trange,prange,yrange=prange,ystyle=1,xstyle=1,/current,/nodata,$ xticklen=.01,xmajor=0,xminor=0,$ ygridstyle=0,yticklen=1,yminor=0,/ylog,$ ytickvalues=[1050,1000,900,800,700,600,500,400,300,200,100],$ ytitle='PRESSURE (hPa)',$ title=title,color='black',font_size=10,position=p_position) ; This writes the xtitle TEMPERATURE t0=text(temp_title_position[0],temp_title_position[1],'TEMPERATURE (!uo!nC)',$ color='black',/normal,alignment=0.5,font_size=10) ; Now need to set up plot space in skewt p1=plot([xmin,xmax],[ymin,ymax],/current,/nodata,xstyle=1,ystyle=1,axis_style=4,$ position=p_position) ;--------------------------------------------------------------------- ; --- DRAW DIAGONAL ISOTHERMS. ;--------------------------------------------------------------------- ; Temp array of -100 to 60 C in 10 deg steps temp=findgen(16)*10.0-100.0 ntemp=n_elements(temp) ; Pressure array from top to bottom p=findgen(prange[0]-prange[1])+prange[1] ; Change y to skewt sy=skewty(p) ; Loop through temperatures so calculate isotherms for i=0,ntemp-1 do begin ; Calculate x in skewt sx=skewtx(temp[i],sy) ; Draw the lines that are within the plot space result=where(sx ge xmin and sx le xmax,count) if count gt 0 then begin l2=polyline(sx[result],sy[result],color='black',target=p1,/data) ; Label at bottom of isotherm xtstring=string(temp[i],format='(I4)') xtstring=strcompress(xtstring,/remove_all) if sy[result[-1]] lt 0 then $ t01=text(sx[result[-1]],sy[result[-1]]-2.0,xtstring,/data,color='black',target=p1,clip=0,alignment=0.5,font_size=10) endif endfor ;--------------------------------------------------------------------- ; --- DRAW DRY ADIABATS. ; --- Iterate in 10 mb increments to compute the x,y points on the ; --- curve. ;--------------------------------------------------------------------- ; Declare adiabat values and pressures where adiabats intersect the ; edge of the skew-T diagram. Refer to a skew-T diagram if necessary. theta=(findgen(24)*10.0)-30.0;temp values from -30 to 170, by 10 degC ntheta=n_elements(theta) ; DRY ADIABATS p=findgen(prange[0]-prange[1])+prange[1] sy=skewty(p) for itheta=0,ntheta-1 do begin dry=tda(theta[itheta], p) sx=skewtx(dry,sy) result=where(sx ge xmin and sx le xmax,count) if count gt 0 then $ l2=polyline(sx[result],sy[result],color='brown',target=p1,/data) endfor ;--------------------------------------------------------------------- ; --- DRAW MOIST ADIABATS UP TO ~ 250hPa ; Declare moist adiabat values and pressures of the tops of the ; moist adiabats. All moist adiabats to be plotted begin at 1050 mb. ;--------------------------------------------------------------------- ; declare pressure range ; convert to plotter coords and declare space for x coords p=reverse(findgen(82)*10.0+240.0);(from = 1050, to = 240, by = -10) npts=n_elements(p) sy=skewty(p) sx=make_array(/double,npts) ; Generating the data for the curves pseudo=[32, 28, 24, 20, 16, 12, 8] NPSEUDO=n_elements(pseudo) ; nrow = npts, ncol = NPSEUDO holdx=make_array(/float,value=0,npts,NPSEUDO) holdy=make_array(/float,value=0,npts,NPSEUDO) for ipseudo=0,NPSEUDO-1 do begin for ilen=0,npts-1 do begin ; satlft is iterative moist=satlft(pseudo[ipseudo], p[ilen]) sx[ilen]=skewtx(moist, sy[ilen]) endfor holdx[*, ipseudo]=sx holdy[*, ipseudo]=sy endfor ; Finally draw the curves. Any curves that extend beyond ; the left axis are clipped. Those curves only get annotated ; at the surface. for ipseudo=0,NPSEUDO-1 do begin ; plot the curves sx=holdx[*, ipseudo] sy=holdy[*, ipseudo] l3=polyline(sx,sy,color='green',target=p1,/data) ; annotate the curves -- at the top moist=satlft(pseudo[ipseudo], 230) labely=skewty(230) labelx=skewtx(moist, labely) if (labelx gt xmin) then begin t3=text(labelx,labely,strcompress(pseudo[ipseudo],/remove_all),color='green',/data,target=p1,font_size=8) ; annotate the curves -- at the surface moist=satlft(pseudo[ipseudo], 1100) labely=skewty(1100) labelx=skewtx(moist, labely) t4=text(labelx,labely,strcompress(pseudo[ipseudo],/remove_all),color='green',/data,target=p1,font_size=8,clip=0,alignment=0.5) endif endfor ;--------------------------------------------------------------------- ; --- DRAW SATURATION MIXING RATIO LINES. ; --- These lines run between 1050 and 400 mb. The 20 line intersects ; --- the sounding below 400 mb, thus a special case is made for it. ; --- The lines are dashed. The temperature where each line crosses ; --- 400 mb is computed in order to get x,y locations of the top of ; --- the lines. ;--------------------------------------------------------------------- mixrat=[20, 12, 8, 5, 3, 2, 1] NMIX=n_elements(mixrat) ; --- Compute y coordinate at the top ; --- (i.e., right end of slanted line) and ; --- the bottom of the lines. ; --- SPECIAL CASE OF MIXING RATIO == 20 yr=skewty(440.) ; y-coord at top (i.e. right) tmix=tmr(mixrat[0], 440.) xr=skewtx(tmix, yr) yl=skewty(1000.) ; y-coord at bot (i.e. left) tmix=tmr(mixrat[0], 1000.) xl=skewtx(tmix, yl) ; dashed line l4=polyline([xl,xr],[yl,yr],color='green',target=p1,/data,linestyle=2) ; We want to stop the mixing ratio lines at 1000 and plot ; the mixing ratio values "in-line" with where the line would continue yl=skewty(1025.) xl=skewtx(tmix, yl) t5=text(xl,yl,strcompress(mixrat[0],/remove_all),color='green',/data,orientation=55,alignment=0.5,target=p1,clip=0,font_size=8) ; --- THE REST OF THE MIXING RATIO marray=make_array(/float,nmix,value=400.0) yr=skewty(marray) tmix=tmr(mixrat,400.) xr=skewtx(tmix,yr) marray=make_array(/float,nmix,value=1000.0) yl=skewty(marray) tmix=tmr(mixrat,1000.) xl=skewtx(tmix, yl) for imix=1,nmix-1 do begin ; dashed line l4=polyline([xl[imix],xr[imix]],[yl[imix],yr[imix]],color='green',target=p1,/data,linestyle=2) endfor marray=make_array(/float,nmix,value=1025.0) yl=skewty(marray) xl=skewtx(tmix, yl) for imix=1,nmix-1 do begin t4=text(xl[imix],yl[imix],strcompress(mixrat[imix],/remove_all),color='green',/data,orientation=55,alignment=0.5,font_size=8,target=p1) endfor ;p0.save,'skewt_axis.png',height=pydim end ;************************************************************************ ; Functions for the skewt_fg ;************************************************************************ ;======================================================================== ; y-coordinate on skew-T, log p diagram given pressure (mb). ; Y-origin at p=1000 mb. ; SKEWTY = 132.182 - 44.061 * ALOG10(PRES) Function skewty,pres skewty=132.18199999999999 - 44.061 * alog10(pres) return,skewty end ;======================================================================== ; Determine x-coordinate on skew-T, log p diagram given ; temperature (C) ; and y-coordinate from FUNCTION SKEWTY. X-origin at T=0c. Function skewtx,temp,ycoord skewtx=0.54000000000000004 * temp + 0.90691999999999995 * ycoord return,skewtx end ;======================================================================== ; This function returns the temperature tda (celsius) on a dry adiabat ; at pressure p (millibars). the dry adiabat is given by potential ; temperature o (celsius). the computation is based on poisson's equation. ; reference stipanuk paper entitled: ; "algorithms for generating a skew-t, log p diagram and computing ; selected meteorological quantities." ; atmospheric sciences laboratory ; u.s. army electronics command ; white sands missile range, new mexico 88002 ; 33 pages, baker, schlatter 17-may-1982 Function tda,o,p ok=o + 273.14999999999998 tdak=ok * ((p * 0.001)^0.28599999999999998) tdak=tdak - 273.14999999999998 return,tdak end ;======================================================================== ; Determine x-coordinate on skew-T, log p diagram given temperature (C) ; and y-coordinate from FUNCTION SKEWTY. X-origin at T=0c. ; ; "algorithms for generating a skew-t, log p ; diagram and computing selected meteorological ; quantities." ; atmospheric sciences laboratory ; u.s. army electronics command ; white sands missile range, new mexico 88002 ; 33 pages, baker, schlatter 17-may-1982 ; this function returns the temperature (celsius) on a mixing ; ratio line w (g/kg) at pressure p (mb). the formula is ; given in table 1 on page 7 of stipanuk (1973). Function tmr,w, p ; initialize constants c1=0.049864645499999999 c2=2.4082965000000001 c3=7.0747499999999999 c4=38.9114 c5=0.091499999999999998 c6=1.2035 x=alog10((w * p)/(622. + w)) tmrk=10^(c1 * x + c2) - c3 + c4 * ((10.^(c5 * x) - c6)^2.) tmrk=tmrk - 273.14999999999998 return,tmrk end ;======================================================================== ; Function satlft,thw,p cta=273.14999999999998 akap=0.28541 pwrp=(p/1000.)^akap tone=(thw + cta) * pwrp - cta eone=wobf(tone) - wobf(thw) rate=1.0 dlt=1.0 while abs(dlt) gt 0.10000000000000001 do begin ttwo=tone - eone * rate pt=(ttwo + cta)/pwrp - cta etwo=pt + wobf(ttwo) - wobf(pt) - thw dlt=etwo * rate rate=(ttwo - tone)/(etwo - eone) tone=ttwo eone=etwo endwhile ttwo=ttwo-dlt return,ttwo end ;======================================================================== ; ; this function calculates the difference of the wet-bulb potential ; temperatures for saturated and dry air given the temperature. ; ;----------------------------------------------------------------------- ; include 'lib_dev:[gudoc]edfvaxbox.for/list' ; baker, schlatter 17-may-1982 original version. ; let wbpts = wet-bulb potential temperature for saturated ; air at temperature t (celsius). let wbptd = wet-bulb potential ; temperature for completely dry air at the same temperature t. ; the wobus function wobf (in degrees celsius) is defined by ; wobf(t) = wbpts-wbptd. ; although wbpts and wbptd are functions of both pressure and ; temperature, their difference is a function of temperature only. ; to understand why, consider a parcel of dry air at tempera- ; ture t and pressure p. the thermodynamic state of the parcel is ; represented by a point on a pseudoadiabatic chart. the wet-bulb ; potential temperature curve (moist adiabat) passing through this ; point is wbpts. now t is the equivalent temperature for another ; parcel saturated at some lower temperature tw, but at the same ; pressure p. to find tw, ascend along the dry adiabat through ; (t,p). at a great height, the dry adiabat and some moist ; adiabat will nearly coincide. descend along this moist adiabat ; back to p. the parcel temperature is now tw. the wet-bulb ; potential temperature curve (moist adiabat) through (tw,p) is wbptd. ; the difference (wbpts-wbptd) is proportional to the heat imparted ; to a parcel saturated at temperature tw if all its water vapor ; were condensed. since the amount of water vapor a parcel can ; hold depends upon temperature alone, (wbptd-wbpts) must depend ; on temperature alone. ; the wobus function is useful for evaluating several thermo- ; dynamic quantities. by definition: ; wobf(t) = wbpts-wbptd. (1) ; if t is at 1000 mb, then t is a potential temperature pt and ; wbpts = pt. thus ; wobf(pt) = pt-wbptd. (2) ; if t is at the condensation level, then t is the condensation ; temperature tc and wbpts is the wet-bulb potential temperature ; wbpt. thus ; wobf(tc) = wbpt-wbptd. (3) ; if wbptd is eliminated from (2) and (3), there results ; wbpt = pt-wobf(pt)+wobf(tc). ; if wbptd is eliminated from (1) and (2), there results ; wbpts = pt-wobf(pt)+wobf(t). ; if t is an equivalent potential temperature ept (implying ; that the air at 1000 mb is completely dry), then wbpts = ept ; and wbptd = wbpt. thus ; wobf(ept) = ept-wbpt. ; this form is the basis for a polynomial approximation to wobf. ; in table 78 on pp.319-322 of the smithsonian meteorological ; tables by roland list (6th revised edition), one finds wet-bulb ; potential temperatures and the corresponding equivalent potential ; temperatures listed together. herman wobus, a mathematician for- ; merly at the navy weather research facility, norfolk, virginia, ; and now retired, computed the coefficients for the polynomial ; approximation from numbers in this table. ; ; notes by t.w. schlatter ; noaa/erl/profs program office ; august 1981 Function wobf,temp x=temp - 20. if x le 0. then begin pol=1. + x * (-0.0088416604999999992 + x * ( $ 0.00014714143000000001 + x * (-9.6719890000000006e-07 + $ x * (-3.2607217000000002e-08 + x * ( $ -3.8598072999999999e-10))))) wbts=15.130000000000001/pol^4 endif else begin pol=1. + x * (0.0036182989000000001 + x * (-1.3603273e-05 + $ x * (4.9618921999999997e-07 + x * ( $ -6.1059364999999998e-09 + x * (3.9401550999999998e-11 + $ x * (-1.2588129e-13 + x * (1.668828e-16))))))) wbts=29.93/pol^4 + 0.95999999999999996 * x - $ 14.800000000000001 endelse return,wbts end