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 * * ---> this version reads the format of file "...*TSPOTINT.txt" * created by IMetOS-II [ 9 Jan 2017, DKM ] * * Compile: * ncargf77 new_skewt.f new_rdsnd2.f -o asc_IMet2_make_skewT * * 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)' selv = 589.0 title = cdate//' UTC at '// + 'Purchase Knob (35.5859N, 83.0729W) sounding (elv=1495 m)' selv = 1495.0 title = cdate//' UTC at '// + 'WW College (35.6091N, 82.4414W) sounding (elv=656 m)' selv = 656.0 title = cdate//' UTC at '// + 'Poga Mntn (36.2526N, 81.9136W) sounding (elv=1018 m)' selv = 1018.0 title = cdate//' UTC at '// + 'UNC Asheville (35.6197N, 82.5675W) sounding (elv=642 m)' selv = 642.0 title = cdate//' UTC at '// + 'UNC Asheville (35.6197N, 82.5675W) sounding (elv=658 m)' selv = 658.0 * * * 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, clin(80)*1 character cdaya(31)*2 * data cnum/'0','1','2','3','4','5','6','7','8','9'/ data cdaya/'01','02','03','04','05','06','07','08','09','10', + '11','12','13','14','15','16','17','18','19','20', + '21','22','23','24','25','26','27','28','29','30', + '31'/ * pi = acos(-1.) * * get date and UTC time from *.SUMMARY file... * * ---> this version reads the format of file "...*SUMMARY.txt" * created by IMetOS-II [ 9 Jan 2017, DKM ] * * * * open(unit=37,file='sounding.STD') open(unit=37,file='sounding.SUM') ifnd=0 10 read(37,110)(clin(j3),j3=1,80) 110 format(80a1) if((clin(1).eq.'L').and.(clin(2).eq.'a').and.(clin(3).eq.'u').and. + (clin(4).eq.'n').and.(clin(5).eq.'c').and.(clin(6).eq.'h'))then igo10=0 ifnd=1 else igo10=1 endif if(igo10.eq.1)goto 10 close(37) * if(ifnd.eq.0)STOP 12009 * jfnd=0 do j3=1,80 if(jfnd.eq.0)then if(clin(j3).eq.':')jfnd=j3 endif enddo * if(jfnd.eq.0)STOP 12019 ind=0 do j3=jfnd,80 if(clin(j3).ne.' ')ind=j3 if(clin(j3).eq.'/')clin(j3)=' ' if(clin(j3).eq.':')clin(j3)=' ' enddo * * year * open(unit=37,file='datntim.dat') write(37,110)(clin(j3),j3=jfnd,ind) close(37) open(unit=37,file='datntim.dat') read(37,*)imo,ida,iyr,ihr,imn,isc close(37) open(unit=37,file='ccyr') write(37,125)iyr 125 format(1x,i4) close(37) open(unit=37,file='ccyr') read(37,127)cyr 127 format(1x,a4) close(37) * * hour (UTC) * * write(*,*)'>>>>',(clin(j4),j4=1,ind),'<<<<' if(clin(ind-1).eq.'P')ihr=ihr+12 imn = imn + nint(float(isc)/60.) * open(unit=37,file='hour.dat') write(37,129)ihr,imn 129 format(1x,i2,1x,i2) close(37) open(unit=37,file='hour.dat') read(37,131)chr(1:2),chr(3:4) 131 format(1x,a2,1x,a2) close(37) if(chr(1:1).eq.' ')chr(1:1)='0' if(chr(3:3).eq.' ')chr(3:3)='0' * * month * if(imo.eq.1)ccmo='01' if(imo.eq.2)ccmo='02' if(imo.eq.3)ccmo='03' if(imo.eq.4)ccmo='04' if(imo.eq.5)ccmo='05' if(imo.eq.6)ccmo='06' if(imo.eq.7)ccmo='07' if(imo.eq.8)ccmo='08' if(imo.eq.9)ccmo='09' if(imo.eq.10)ccmo='10' if(imo.eq.11)ccmo='11' if(imo.eq.12)ccmo='12' * * day * cday = cdaya(ida) * cdat=cyr//' '//ccmo//' '//cday//' '//chr(1:2)//':'//chr(3:4) * write(*,*)' ' write(*,*)' sounding date=',cdat,'<<<<' * iredo=0 17 iskip=0 ibop=0 ilvl=0 do ihdr=1,3 read(iunit,*) enddo 20 ilvl=ilvl+1 tdc = -9999.0 read(iunit,*,end=500)time,pmb,tc,rh,ws,wd,zm selv = 658.0 * ASL zm = zm + selv * 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 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 xs(1,6,ilvl)=wd * wd in knots originally [convert to m/s] * xs(1,7,ilvl)=ws*0.514444 xs(1,7,ilvl)=ws * 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