program process_snd c c program to process hires sounding data obtained from ASPEN in c EOL or CLS format c creates a uniform pressure dataset ... upapi_stnid c c gfortran create_prs_L4.f -o create_prs_L4.x c c 1) read in data c 2) interpolate data to specific uniform vertical pressure levels c 3) put qc flags on data (all flags are set as good) real xmis, zmis, ixmis parameter (ip=1025, xmis=-999., zmis=-999., ixmis=-999.) real p(ip), z(ip), tc(ip), td(ip), uc(ip), vc(ip) real fuc(ip), fvc(ip), xlni(ip), xlti(ip) real psfc(1000) integer qp(ip), qh(ip), qt(ip), qd(ip), qu(ip) integer istnid, nf parameter (m=10000) real xlon, xlat, alt integer iyr, mon, idy, ihr, min integer imm, idd, ihh common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 2 qpp(m), qtt(m), qhh(m), quu(m), qvv(m), xln(m), xlt(m) character fname*33, cid*5, fout*12, resp*1 integer igood, iques, idoub, ibad, iint, iest, imis, iunc data igood, iques, ibad, iint, iest, iunc, imis 2 / 1, 2, 4, 6, 7, 8, 9/ c ************************************************************ c specify constants cid = ' ' print *, 'read in stnid (e.g., 46810) ...>' read(5,*) istnid print *, 'istnid = ', istnid print *, 'Is this correct (y/n) ...' read(5,*) resp if(resp .eq. 'n' .or. resp .eq. 'N') then print *, 'incorrect cid and istnid' stop endif print *, 'if humidity has been corrected type 1, else 0 ...>' read(5,*) ihumcor write(fout, 8) istnid 8 format('upapi_',i5) open (20, file=fout, status='unknown', form='formatted') c c specify resolution of dataset in mb c print *, 'Specify desired uniform resolution in mb ' print *, '(1 mb or greater ...)' read (5,*) pinc r = 287.054 g = 9.806 pi = acos(-1.) degran = pi / 180. c contains list of files to be processed open(2, file='files_hires', form='formatted', status='old') nf = 0 5 continue c set quality control flags to good at all levels do l=1,ip qp(l) = igood qh(l) = igood qt(l) = igood qd(l) = igood qu(l) = igood enddo c read in file read(2,'(a)',end=999) fname nf = nf + 1 print *, 'processing data for ', nf, fname c open data files c file 8 contain input hi-res data open(8, file=fname, form='formatted', status='old') c c read in one data for one time from file 8 c c select format(cls or eol) that input data is in by c uncommenting appropriate line c call readcls(8, nl) call readeol(8, nl) c c get pressures to interpolate the data to c psfc(nf) = p0(1) p(1) = p0(1) p(2) = ifix(p0(1)/pinc) * pinc if(p(1) .eq. p(2)) then print 233, mon, idy, ihr 233 format(40x, 'check ', 3i3) p(2) = p(1) - pinc endif n = 2 60 n = n + 1 p(n) = p(n-1) - pinc if(p(n) .gt. p0(nl)) go to 60 np = n - 1 c get data at interpolated levels call interp(p0, z0, nl, p, z, np, zmis) call interp(p0, t0, nl, p, tc, np, xmis) call interp(p0, d0, nl, p, td, np, xmis) call layave(p0, pinc, u0, nl, p, uc, np, xmis) call layave(p0, pinc, v0, nl, p, vc, np, xmis) call fil121(uc, np, fuc, ixmis) call fil121(vc, np, fvc, ixmis) c call interp(p0, u0, nl, p, uc, np, xmis) c call interp(p0, v0, nl, p, vc, np, xmis) call interp(p0, xln, nl, p, xlni, np, xmis) call interp(p0, xlt, nl, p, xlti, np, xmis) c compute date/time tag call dattim(iyr, mon, idy, ihr, min, imm, idd, ihh) write(11,222) fname, mon, idy, ihr*100+min 222 format(a19,3i5) c print out header information write(20,700) 700 format(' STN DATE GMT HTS LAT LONG ID') c write(20,701) cid, iyr,imm,idd,ihh*100, nint(alt), xlat, c 2 xlon, istnid write(20,701) cid, iyr,mon,idy,ihr*100+min, nint(alt), xlat, 2 xlon, istnid 701 format(a4, 4x, 3i2.2, 1x, i4.4, 1x, i5, 2(1x,f8.2), i7) if (ihumcor .eq. 0) then write(20,702) np else write(20,704) np endif 702 format('NLVL =',i4,' COMMENTS:') 704 format('NLVL =',i4,' COMMENTS: humidity corrected') write(20,703) 703 format(' P HT TC TD DIR SPD QP QH', 2 ' QT QD QW LON LAT') do it=1,np if(p(it) .eq. xmis) qp(it) = imis if(z(it) .eq. xmis) qh(it) = imis if(tc(it) .eq. xmis) qt(it) = imis if(td(it) .eq. xmis) qd(it) = imis if(fuc(it) .ne. xmis .and. fvc(it) .ne. xmis) then call spd_dir(fuc(it), fvc(it), spd, dir, ixmis) qu(it) = igood else spd = ixmis dir = ixmis qu(it) = imis endif if(xlni(it) .eq. ixmis) then xlni(it) = xlni(it-1) endif if(xlti(it) .eq. ixmis) then xlti(it) = xlti(it-1) endif write (20,200) p(it),z(it),tc(it),td(it), dir, spd, 2 qp(it), qh(it), qt(it), qd(it), qu(it), xlni(it), xlti(it) enddo 200 format(6f8.1,1x,5i3,2f8.2) go to 5 999 continue ntot = nf - 1 print *, 'number of sondes processed ', ntot call mean (psfc, ntot, apsfc, xmis, nn) call PFMIMA_mis( ntot, psfc, pMIN, pMAX, xmis ) print *, 'mean surface pressure = ', apsfc print *, 'number of obs = ', nn print *, 'pmin, pmax = ', pmin, pmax end subroutine fil121(fld, n, ffld, xmis) c c filter data with 1-2-1 filter c real fld(n), ffld(n), xmis c ******************************************************************* ffld(1) = fld(1) ffld(n) = fld(n) do i=2,n-1 ip = i + 1 im = i - 1 if(fld(ip).ne.xmis .and. fld(i).ne.xmis 2 .and. fld(im).ne.xmis) then ffld(i) = 0.25*(fld(im)+fld(ip)) + 0.5*fld(i) else ffld(i) = fld(i) endif enddo return end subroutine dattim(iyr, mon, idy, ihr, min, imm, idd, ihh) integer iyr, mon, idy, ihr, min, imm, idd, ihh c add one hour to sonde time (sondes are typically launched 1 hr c prior to time that obs are needed) ihh = ihr + nint(float(min)/60.) + 1 imm = mon idd = idy if(ihh .eq. 24) then ihh = 1 idd = idy + 1 if(mon .eq. 6 .and. idd .gt. 30) then imm = 7 idd = 1 elseif(mon .eq. 7 .and. idd .gt. 31) then imm = 8 idd = 1 elseif(mon .eq. 8 .and. idd .gt. 31) then imm = 9 idd = 1 endif elseif(ihh .gt. 25) then print *, 'error in dattim', ihh, imm, idd stop endif return end subroutine spd_dir(u, v, spd, dir, xmis) c c convert u and v wind components into speed and direction c pi = acos(-1.) ran2deg = 180. / pi if(u.ne.xmis .and. v.ne.xmis) then spd = sqrt (u*u + v*v) if (spd .eq. 0.0) then dir = 0. else if(u .eq. 0.) then if (v .ge. 0.) then dir = 180. else dir = 360. endif elseif(v .eq. 0.) then if (u .ge. 0.) then dir = 270. else dir = 90. endif else dir = acos(-v/spd) * ran2deg if(u .gt. 0) dir = 360. - dir endif endif else spd = xmis dir = xmis endif return end subroutine layave(p, pinc, fld, nl, pi, fldi, np, xmis) c this routine computes a smooth profile by taking layer averages c c maximum interval over which to do linear interpolation c real pmax, imis parameter (pmax = 50., imis=-999.) c raw fields real p(nl), fld(nl) c interpolated fields real pi(np), fldi(np) c missing value flag real xmis c ********************************************************************** c set surface value to first raw value if(fld(1) .ne. xmis) then fldi(1) = fld(1) else fldi(1) = imis endif do 70 il=2,np c set field to missing but fill in if possible fldi(il) = imis pa = pi(il) - pinc/2. if(pa .lt. p(nl)) then pa = p(nl) endif pb = pi(il) + pinc/2. if(pb .gt. p(1)) then pb = p(1) endif sum = 0. nobs = 0 do iraw = 1,nl if(p(iraw) .lt. pb .and. p(iraw).gt.pa) then if(fld(iraw) .ne. xmis) then sum = sum + fld(iraw) nobs = nobs + 1 endif endif enddo if(nobs .gt. 0) then fldi(il) = sum / float(nobs) endif 70 continue return end subroutine interp(p, fld, nl, pi, fldi, np, xmis) c this routine linearly interpolates a field in log pressure c c maximum interval over which to do linear interpolation c real pmax, imis parameter (pmax = 50., imis=-999.) c raw fields real p(nl), fld(nl) c interpolated fields real pi(np), fldi(np) c missing value flag real xmis c ********************************************************************** c set surface value to first raw value fldi(1) = fld(1) do 70 il=2,np c set field to missing but fill in if possible fldi(il) = imis c this loop gets fld if interpolated pressure level is a raw data level do iraw = 1,nl if(pi(il) .eq. p(iraw)) then if(fld(iraw) .eq. xmis) go to 10 fldi(il) = fld(iraw) c print *, p(il), fldi(il) go to 70 endif enddo c this loop get fld if interpolated pressure level is not a raw data level c then interpolated value in linearly interpolated in log of pressure c find iplv - pressure level in raw data closest to p(il) 10 continue pdif = 1.e4 do iraw = 1,nl dp = abs(p(iraw) - pi(il)) c if(iraw .lt. 20 )print *, iraw, p(iraw), pi(il), dp, pdif, iplv if(dp .lt. pdif) then pdif = dp iplv = iraw endif enddo c print *, il, pi(il), p(iplv) c find first level below pi(il) in which data is not missing do 20 ip=iplv,1,-1 if(p(ip) .lt. pi(il)) go to 20 if(fld(ip) .ne. xmis) then ipb = ip go to 30 endif 20 continue go to 70 30 continue c find first level above pi(il) in which data is not missing do 40 ip=iplv,nl if(p(ip) .gt. pi(il)) go to 40 if(fld(ip) .ne. xmis) then ipa = ip go to 50 endif 40 continue go to 70 50 continue if(abs(p(ipb)-p(ipa)) .gt. pmax) go to 70 ap = alog(pi(il)) apb = alog(p(ipb)) apa = alog(p(ipa)) c print *, il, iplv, ipb, ipa c print *, pi(il), p(ipb), p(ipa) c print *, fld(ipb), fld(ipa) sl = (ap - apb) / (apa - apb) fldi(il) = fld(ipb) + sl*(fld(ipa) - fld(ipb)) c print *, fld(ipb), fldi(il) 70 continue return end subroutine readeol(iun, lvl) c c read in data for one time c assuming data is in NCAR CLass EOL parameter (m=10000, pmis=-999.) real time, rh, wc, wspd, wdir, dz, xlnt, xltt integer iyr, mon, idy, ihr, min, ih, im real sec, gpsalt common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 2 qp(m), qt(m), qh(m), qu(m), qv(m), xln(m), xlt(m) character hline*100, dline*135 c ************************************************************* c c routine to read in uncorrected souding data c iflag = 0 10 read (iun,'(a100)',end=100) hline if (hline(1:4) .eq. 'Data') then c skip over a few lines read(iun, '(a100)',end=100) hline read(iun, '(a100)',end=100) hline read(iun, '(a100)',end=100) hline c read in lon, lat, alt info read(iun, '(a100)',end=100) hline read (hline(56:65),*) xlon read (hline(79:87),*) xlat read (hline(91:97),*) alt c read in date / time information read (iun,'(a100)') hline read (hline(46:47),*) iyr read (hline(50:51),*) mon read (hline(54:55),*) idy read (hline(58:59),*) ihr read (hline(61:62),*) min c read remainder of header lines do i=1,8 read (iun,'(a82)') hline enddo else go to 10 endif C------- read in data --------------------------------------------------- lvl = 0 nl = 0 20 read (iun,'(a133)',end=100) dline nl = nl + 1 read (dline,*) time, ih, im, sec, pres if(pres .ne. pmis) then lvl = lvl + 1 if(lvl .gt. m) go to 100 read (dline,*) time, ih, im, sec, p0(lvl), t0(lvl), d0(lvl), 2 rh, u0(lvl), v0(lvl), wspd, wdir, dz, z0(lvl), xlnt, 3 xltt, gpsalt if(xlnt .eq. -999.) then if(lvl .eq. 1) then xln(lvl) = xlon else xln(lvl) = xln(lvl-1) endif else xln(lvl) = xlnt endif if(xlt(lvl) .eq. -999.) then if(lvl .eq. 1) then xlt(lvl) = xlat else xlt(lvl) = xlt(lvl-1) endif else xlt(lvl) = xltt endif c print 702, p0(lvl), z0(lvl), t0(lvl), d0(lvl), c 2 u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(lvl) c print 702, p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl) c write(49,702) p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl) endif go to 20 100 continue print 222, xlon, xlat, alt, 2 iyr, mon, idy, ihr, min,nl,lvl 222 format (2f7.2,f4.1,7i5,f7.2) 702 format (9(1x,f7.1)) return end subroutine readcls(iun, lvl) c c read in data for one time c assuming data is in NCAR CLass format parameter (m=10000, pmis=9999.) real time, rh, wc, wspd, wdir, dz, rng, xlnt, xltt integer iyr, mon, idy, ihr, min common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 2 qp(m), qt(m), qh(m), qu(m), qv(m), xln(m), xlt(m) character hline*100, dline*100 c ************************************************************* c c routine to read in uncorrected souding data c iflag = 0 10 read (iun,'(a100)',end=100) hline if (hline(1:4) .eq. 'Data') then c skip over a few lines read(iun, '(a100)',end=100) hline read(iun, '(a100)',end=100) hline c read in lon, lat, alt info read(iun, '(a100)',end=100) hline c read for 46734 read (hline(61:69),*) xlon read (hline(71:78),*) xlat read (hline(80:82),*) alt c read (hline(63:73),*) xlon c read (hline(74:82),*) xlat c read (hline(84:85),*) alt c read in date / time information read (iun,'(a100)') hline read (hline(38:39),*) iyr read (hline(42:43),*) mon read (hline(46:47),*) idy read (hline(50:51),*) ihr read (hline(53:54),*) min c read remainder of header lines do i=1,10 read (iun,'(a82)') hline enddo else go to 10 endif C------- read in data --------------------------------------------------- lvl = 0 nl = 0 20 read (iun,'(a100)',end=100) dline nl = nl + 1 read (dline,*) time, pres if(pres .ne. pmis) then lvl = lvl + 1 if(lvl .gt. m) go to 100 read (dline,*) time, p0(lvl), t0(lvl), d0(lvl), rh, 2 u0(lvl), v0(lvl), wspd, wdir, dz, xlnt, 3 xltt, rng, az, z0(lvl) if(xlnt .ge. 999.) then if(lvl .eq. 1) then xln(lvl) = xlon else xln(lvl) = xln(lvl-1) endif else xln(lvl) = xlnt endif if(xlt(lvl) .ge. 999.) then if(lvl .eq. 1) then xlt(lvl) = xlat else xlt(lvl) = xlt(lvl-1) endif else xlt(lvl) = xltt endif c print 702, p0(lvl), z0(lvl), t0(lvl), d0(lvl), c 2 u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(lvl) c print 702, p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl) c write(49,702) p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl) endif go to 20 100 continue print 222, xlon, xlat, alt, 2 iyr, mon, idy, ihr, min,nl,lvl,p0(lvl) 222 format (2f7.2,f4.1,7i5,f7.2) 702 format (9(1x,f7.1)) return end SUBROUTINE MEAN(DATA, N, EZMEAN, XMIS, nn) C C THIS CALCULATES THE MEAN OF A SET OF DATA C REAL DATA(1) SUM = 0. nn = 0 DO 10 I=1,N if(data(i) .ne. xmis) then SUM = SUM + DATA(I) nn = nn + 1 endif 10 CONTINUE if(nn .gt. 0) then EZMEAN = SUM/FLOAT(NN) else EZMEAN = xmis endif RETURN END SUBROUTINE PFMIMA_mis( N, X, XMIN, XMAX, xmis ) INTEGER N REAL X(N), XMIN, XMAX C C Finds the minimum and maximum values in the vector X of length N C INTEGER I REAL BIG PARAMETER ( BIG = 1.0E29 ) C XMIN = +BIG XMAX = -BIG DO 10 I=1,N if(x(i) .eq. xmis) go to 10 XMIN = MIN( XMIN, X(I) ) XMAX = MAX( XMAX, X(I) ) 10 CONTINUE RETURN END