program assign_qc c c program to process hires sounding data through a series c of objective qc checks where data is in uniform height resolution. c c gfortran assign_qcf_hgt_L4.f -o assign_qcf_hgt_L4.x c c 1) read in data c 2) read in stats created by create_qcstats_prs_L4.x c 3) assign QC flags real xmis parameter (maxz=10000, ntmx=250, xmis=-999.) real pres(maxz), zs(ntmx) real p(maxz), z(maxz), tc(maxz), td(maxz), uc, vc real dir(maxz), spd(maxz), xln(maxz), xlt(maxz) integer qp(maxz), qh(maxz), qt(maxz), qd(maxz), qw(maxz) integer nt, nz, nu, nv, nd integer idate(ntmx), itime(ntmx) real pave(maxz), psd(maxz) real tave(maxz), tsd(maxz) real tdave(maxz), tdsd(maxz) real uave(maxz), usd(maxz) real vave(maxz), vsd(maxz) real xlon, xlat, alt integer iyr, mon, idy, ihr, min integer imm, idd, ihh character fin*11, fstat*15, fout*11, header(5)*64, resp*1 character ferr*13 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 fac1 - number of standards to allow before questionable flag is set fac1 = 3. c fac2 - number of standards to allow before bad flag is set fac2 = 6. npq = 0 npb = 0 ntq = 0 ntb = 0 ndq = 0 ndb = 0 nwq = 0 nhb = 0 ntdb = 0 ntdq = 0 c specify constants print *, 'read in wmo number of site to process ...>' read(5,'(i5)') istnid print *, 'specify zinc of data ' read(5, *) zinc print *, 'enter stn elev ...>' read(5,*) zstart c print *, ' istnid = ', istnid c print *, 'Is this correct (y/n) ...' c read(5,*) resp c if(resp .eq. 'n' .or. resp .eq. 'N') then c print *, 'incorrect istnid' c stop c endif c initialize fields write(fin, 8) istnid 8 format('upazi_',i5) write(fout, 18) istnid 18 format('upaqf_',i5) write(ferr, 14) istnid 14 format('qchecks_',i5) write(fstat, 9) istnid 9 format('stats_hgt_',i5) open (11, file=ferr, status='unknown', form='formatted') open (20, file=fin, status='old', form='formatted') open (22, file=fstat, status='old', form='formatted') open (24, file=fout, status='unknown', form='formatted') k = 0 3 k = k + 1 read(22,222,end=4) kdum, z(k), pave(k), psd(k), nz, tave(k), 2 tsd(k), ntc,tdave(k), tdsd(k), nd, uave(k), usd(k), 3 nu, vave(k), vsd(k), nv go to 3 4 max = k - 1 print *, 'read in ', max, ' lines of stats ' c read header information nt = 0 5 nt = nt + 1 do il=1,4 read(20,'(a64)',end=50) header(il) write(24,'(a64)') header(il) enddo read(header(2)(9:14),*) idate(nt) read(header(2)(16:19),*) itime(nt) read(header(3)(7:10),'(i4)') nlv write(11,*) 'checks for: ', idate(nt), itime(nt), nlv c print *, idate(nt), itime(nt), nlv do k=1,nlv read (20,200,end=50) p(k), z(k), tc(k), td(k), dir(k), spd(k), 2 qp(k), qh(k), qt(k), qd(k), qw(k), xln(k), xlt(k) if(k .eq. 1) then izindex = 1 else izindex = (z(k) - zstart)/zinc + 2 endif c print *, z(k), izindex, p(k), pave(izindex), psd(izindex) c check pressure dif = abs(p(k) - pave(izindex)) if(p(k) .eq. xmis) then qp(k) = imis elseif(dif .gt. fac2*psd(izindex)) then qp(k) = ibad npb = npb + 1 write(11,*) 'z bad', p(k), z(k), pave(izindex), dif elseif(dif .gt. fac1*psd(izindex)) then qp(k) = iques npq = npq + 1 write(11,*) 'z ques', p(k), z(k), pave(izindex), dif endif c check temperature dif = abs(tc(k) - tave(izindex)) if(tc(k) .eq. xmis) then qt(k) = imis elseif(dif .gt. fac2*tsd(izindex)) then write(11,*) 't bad', p(k), tc(k), tave(izindex), dif qt(k) = ibad ntb = ntb + 1 elseif(dif .gt. fac1*tsd(izindex)) then write(11,*) 't ques', p(k), tc(k), tave(izindex), dif qt(k) = iques ntq = ntq + 1 endif c check dew point c double dew point check criteria dif = abs(td(k) - tdave(izindex)) if(td(k) .eq. xmis) then qd(k) = imis elseif(dif .gt. 2*fac2*tdsd(izindex) .or. td(k).gt.tc(k)) then write(11,*) 'td bad', p(k), td(k), tdave(izindex), dif qd(k) = ibad ndb = ndb + 1 elseif(dif .gt. 2*fac1*tdsd(izindex)) then write(11,*) 'td ques', p(k), td(k), tdave(izindex), dif qd(k) = iques ndq = ndq + 1 endif c check winds ******************************************************* call spdr_uv (spd(k), dir(k), uc, vc, xmis) difu = abs(uc - uave(izindex)) difv = abs(vc - vave(izindex)) if(spd(k) .eq. xmis .or. dir(k) .eq. xmis) then qw(k) = imis elseif(difu .gt. fac2*usd(izindex) .or. 2 difv .gt. fac2*vsd(izindex) ) then qw(k) = ibad nwb = nwb + 1 write(11,*) 'winds bad', p(k), uc, vc, difu, difv elseif(difu .gt. fac1*usd(izindex) .or. 2 difv .gt. fac1*vsd(izindex) ) then qw(k) = iques nwq = nwq + 1 write(11,*) 'winds ques', p(k), uc, vc, difu, difv endif enddo c c vertical consistency checks c c on geopotential height c nhs = 0 do 70 l=1,nlv-1 if(z(l+1) .eq. xmis .or. z(l).eq.xmis) go to 70 if(z(l+1) .lt. z(l)) then write(11,130) p(l), z(l+1), z(l) nhs = nhs + 1 endif 70 continue c c on pressure c nmp = 0 do 80 l=1,nlv-1 if(p(l+1) .eq. xmis .or. p(l).eq.xmis) go to 80 if(p(l+1) .gt. p(l)) then write(11,135) p(l) nmp = nmp + 1 endif 80 continue c c on lapse rate c nlr = 0 do 90 l=1,nlv-1 if(tc(l+1).ne.xmis .and. tc(l).ne.xmis .and. 2 z(l+1).ne.xmis .and. z(l).ne.xmis) then dtdz = (tc(l+1) - tc(l)) / (z(l+1)-z(l)) * 1000. dp = p(l) - p(l-1) if(l.eq.1 .or. l.eq.2) then if(abs(dtdz) .gt. 50.) then write(11,140) p(l), dtdz, tc(l), tc(l+1), p(l), p(l+1) nlr = nlr + 1 endif go to 90 endif if(dp.le.3 .and. abs(dtdz) .le. 30.) then go to 90 elseif(dtdz .lt. -15.) then write(11,140) p(l), dtdz, tc(l), tc(l+1), p(l), p(l+1) nlr = nlr + 1 elseif(dtdz .gt. 15. .and. p(l).lt.900. 2 .and. p(l).gt.100.) then write(11,140) p(l), dtdz, tc(l), tc(l+1), p(l), p(l+1) nlr = nlr + 1 elseif(dtdz .gt. 30.) then write(11,140) p(l), dtdz, tc(l), tc(l+1), p(l), p(l+1) nlr = nlr + 1 endif else go to 90 endif 90 continue c on sfc dew point if(td(2) .ne. xmis .and. td(1).ne.xmis) then dtd = td(1) - td(2) else dtd = xmis endif if(dtd .gt. 6) then ntdb = ntdb + 1 qd(1) = ibad print *, idate(nt), itime(nt), dtd, 'bad' else if(dtd .gt. 3) then ntdq = ntdq + 1 qd(1) = iques print *, idate(nt), itime(nt), dtd, 'ques' endif do k=1,nlv write (24,200) p(k), z(k), tc(k), td(k), dir(k), spd(k), 2 qp(k), qh(k), qt(k), qd(k), qw(k), xln(k), xlt(k) enddo c stop go to 5 50 continue 200 format(6f8.1,1x,5i3,2f8.2) 999 continue print *, 'processed ', nt -1, 'snds' print *, 'h: ques, bad ', npq, npb print *, 't: ques, bad ', ntq, ntb print *, 'td: ques, bad ', ndq, ndb print *, 'w: ques, bad ', nwq, nwb print *, 'hydrostatic check ', nhs print *, 'mono pres check ', nmp print *, 'lapse rate check ', nlr print *, 'sfc td check ', n 130 format(' ***** z not increasing at ',f6.0,' value ... ',2f7.2) 135 format(' ***** p not decreasing at ',f6.0) 140 format('check lapse rate at ',4f7.2,2f7.0) 222 format(i5,f7.0,1x,3(f6.1,f6.2,i4),2(f7.1,f6.1,i4)) 223 format(3i6,f8.1) end SUBROUTINE VARC(DATA, N, VAR, SDEV, XMEAN, XMIS) c c routine to compute the standard deviation and variance c REAL DATA(1), VAR, XMEAN M = 0 SUM = 0. VAR = -99. SDEV = -99. DO 20 I=1,N IF(DATA(I) .EQ. XMIS) GO TO 20 M = M + 1 SUM = SUM + (DATA(I) - XMEAN)**2 20 CONTINUE if(m .gt. 0) then VAR = SUM / FLOAT(M-1) SDEV = SQRT(VAR) endif RETURN END SUBROUTINE MEAN(DATA, N, EZMEAN, XMIS, nn) C C THIS CALCULATES THE MEAN OF A SET OF DATA C REAL DATA(1) c ********************************************* SUM = 0. nn = 0 ezmean = xmis 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) endif RETURN 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 spdr_uv (spd, dir, u, v, xmis) c c convert u and v wind components into speed and direction c pi = acos(-1.0) degran = pi / 180. u = xmis v = xmis if (dir.gt.360. .or. dir.lt.0. .or. spd.lt.0.) then return else u = -1.* spd * sin(dir*degran) v = -1.* spd * cos(dir*degran) 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 readupa(iun, lvl) c c read in data for one time parameter (m=10000, pmis=9999.) real time, rh, wc, wspd, wdir, dz, xln, xlt, rng 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) character hline*130, dline*130 c ************************************************************* c c routine to read in uncorrected souding data c iflag = 0 10 read (iun,'(a82)',end=100) hline if (hline(1:10) .eq. 'Launch Loc') then c read in lon, lat, alt info read (hline(61:70),*) xlon read (hline(72:78),*) xlat read (hline(81:81),*) alt c read in date / time information read (iun,'(a82)') 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,'(a130)',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), wc, wspd, wdir, dz, xln, xlt, rng, 3 z0(lvl), qp(lvl), qt(lvl), qh(lvl), qu(lvl), qv(lvl) c print 702, p0(lvl), z0(lvl), t0(lvl), d0(lvl), c 2 u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(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