subroutine rdsnd(iunit,pres,temp,dwpt,spd,dir,nn,kl,title,badval) * * This subprogam serves three purposes (porpii ?) * * [1] read in NWFS actual soundings * [2] convert old RH/Td to corrected RH/Td * [3] pass corrected sounding to skew-T plotting routine * * parameter(niops=1,nsta=1,nsndg=1,npar=7,nlvl=20000) parameter(nens=9) * real pres(nn), temp(nn), dwpt(nn), spd(nn), dir(nn) real xs(nsta,npar,nlvl), xn(nsta,npar,nlvl) * integer nsnd(niops), klvl(nsndg), mlvl(nsndg) * character cdat(nsndg)*16 character cdate*16, title*80 * data cdat/1*'8888888888888888'/ * if(nlvl.ne.nn)then write(*,*)' nn from MAIN .ne. nlvl from RDSND' write(*,*)' nn=',nn write(*,*)' nlvl=',nlvl STOP 30232 endif * badval = -9999. * do iiop=1,1 ! loop over number of IOPs * call rdsonde(iunit,cdate,ilvl,nsta,npar,nlvl,xs) * title = cdate//' UTC at '// + 'ETS Univ (36.2944N, 82.3706W) sounding (elv=589 m)' title = cdate//' UTC at '// + 'Purchase Knob (35.5859N, 83.0729W) sounding (elv=1495 m)' title = cdate//' UTC at '// + 'WW College (35.6091N, 82.4414W) sounding (elv=656 m)' title = cdate//' UTC at '// + 'Poga Mntn (36.2526N, 81.9136W) sounding (elv=1018 m)' title = cdate//' UTC at '// + 'UNC Asheville (35.6197N, 82.5675W) sounding (elv=642 m)' title = cdate//' UTC at '// + 'UNC Asheville (35.6197N, 82.5675W) sounding (elv=658 m)' * * * write(*,*)' converting RH to dewpoint, ilvl=',ilvl call RHfix(ilvl,nsta,npar,nlvl,xs,xn) * do kk=1,ilvl pres(kk) = xs(1,2,kk) temp(kk) = xs(1,3,kk) dwpt(kk) = xn(1,5,kk) spd(kk) = xs(1,7,kk) dir(kk) = xs(1,6,kk) enddo * kl = ilvl-1 * write(*,*)' ' write(*,*)' ' write(*,*)' returning CONVERTED sounding file, ilvl=',ilvl * enddo ! loop over number of IOPs * return end subroutine rdsonde(iunit,cdat,ilvl,nsta,npar,nlvl,xs) * * iMet 3050 *.LOG file format... * * c1 = T (degC), c2 = Tv (K), c3 = RH (%), c4 = Pres (hPa), c5 = GeoHt (m), * c6 = RDF elev angle (rad), c7 = RDF azim angle (rad), c8 = Wind Dir (rad), * c9 = Wind Speed (knots), c10 = slant range (m), c11 = status flag, * c12 = lat (deg), c13 = long (deg), c14 = Traw (degC), * c15 = geometric height (m MSL), c16 = record time (s), c17 = height status * real xs(nsta,npar,nlvl) * character cdat*16, cyr*4, ccmo*2, cday*2, chr*4, cmo*2 * character inline*120, cline*80, cnum(10)*1 * data cnum/'0','1','2','3','4','5','6','7','8','9'/ * pi = acos(-1.) * * get date and UTC time from *.STD file... * open(unit=37,file='sounding.STD') read(37,110)cday,cmo,cyr,chr(1:2),chr(3:4) 110 format(31x,a2,1x,a2,1x,a4,8x,a2,1x,a2) * read station elevation in STD ??? close(37) * ccmo=cmo if(cmo.eq.' 1')ccmo='01' if(cmo.eq.' 2')ccmo='02' if(cmo.eq.' 3')ccmo='03' if(cmo.eq.' 4')ccmo='04' if(cmo.eq.' 5')ccmo='05' if(cmo.eq.' 6')ccmo='06' if(cmo.eq.' 7')ccmo='07' if(cmo.eq.' 8')ccmo='08' if(cmo.eq.' 9')ccmo='09' cdat=cyr//' '//ccmo//' '//cday//' '//chr(1:2)//':'//chr(3:4) * write(*,*)' ' write(*,*)' sounding date=',cdat,'<<<<' * iredo=0 17 iskip=0 ibop=0 ilvl=0 20 ilvl=ilvl+1 tdc = -9999.0 read(iunit,*,end=500)tc,tv,rh,pmb,zm,rea,raa,wd,ws * read(iunit,*,end=500)xmin,xsec,arat,zm,pmb,tc,rh,tdc,wd,ws * read(iunit,250,end=500)cline 250 format(a80) if(ilvl.gt.nlvl)then write(*,*)' ' write(*,*)' ilvl.GT.nlvl',ilvl,nlvl iredo=1 ibop=1 endif * if(ibop.eq.1)then close(iunit) open(iunit,file='sounding.txt') goto 17 endif * if((iredo.eq.1).and.(pmb.lt.500.))then if(iskip.eq.0)then ilvl=ilvl-1 iskip=1 else iskip=0 endif endif if(iskip.eq.1)goto 20 * xs(1,1,ilvl)=zm xs(1,2,ilvl)=pmb xs(1,3,ilvl)=tc xs(1,4,ilvl)=rh xs(1,5,ilvl)=tdc * wd in radians originally [convert to degrees] xs(1,6,ilvl)=wd*180./pi * wd in knots originally [convert to m/s] xs(1,7,ilvl)=ws*0.514444 * skew-T plot has 100 mb level as its limit if(pmb.lt.100.)goto 500 * goto 20 * 500 ilvl=ilvl-1 * 500 continue * return end subroutine RHfix(ilvl,nsta,npar,nlvl,xs,xn) ! convert from RH to dewpoint temp for iMet sondes * real xs(nsta,npar,nlvl), xn(nsta,npar,nlvl) * * data aa/-8.3536441544146420E-006/ * data bb/8.8877859894645788E-003/ * data cc/1.000000000000000/ * do kk=1,ilvl pmb = xs(1,2,kk) tc = xs(1,3,kk) rhold = xs(1,4,kk) tdold = xs(1,5,kk) rhnew=rhold ** tdnew=tdold * * B = (ln(RH / 100) + ((17.27 * T) / (237.3 + T))) / 17.27 * * D = (237.3 * B) / (1 - B) * * where: * T = Air Temperature (Dry Bulb) in Centigrade (C) degrees * RH = Relative Humidity in percent (%) * B = intermediate value (no units) * D = Dewpoint in Centigrade (C) degrees * * http://ag.arizona.edu/azmet/dewpoint.html * B = (alog(rhnew/100) + ((17.27 * tc) / (237.3 + tc))) / 17.27 tdnew = (237.3 * B) / (1.0 - B) * do ipar=1,7 xn(1,ipar,kk)=xs(1,ipar,kk) enddo xn(1,4,kk)=rhnew xn(1,5,kk)=tdnew enddo * return end