pro edit_EDT,dd=dd,mm=mm,yy=yy,hh=hh ; +---------------------------+ ; | | ; | AUTHOR: Henry Selkirk | ; | DATE: 25 Aug 2005 | ; | | ; +---------------------------+ ; First create input filename from arguments passed from calling program cdd = string(dd) if dd lt 10 then cdd = '0' + cdd cmm = string(mm) if mm lt 10 then cmm = '0' + cmm cyy = string(yy) if yy lt 10 then cyy = '0' + cyy chh = string(hh) if hh lt 10 then chh = '0' + chh txtfile = string(cdd) + string(cmm) + string(cyy) + string(chh) + 'z.txt' txtfile = strcompress(txtfile,/remove_all) print,txtfile ; Next assemble the path - Here we assume that you are running from a directory in which the ; subdirectories 'data' and 'EDT_archive' reside yyyyPathmm = strcompress('20' + cyy + '-' + cmm,/remove_all) print,yyyyPathmm hhZ = strcompress( chh + 'Z',/remove_all) inPath = 'data' + '/' + yyyyPathmm + '/' + hhZ + '/' + txtfile print,inPath ; Open the input file openr,12,inPath ; ************************************************************************* ; The message data are written into the special comment lines which come BEORE the normal header data that ; are written into the normal comment lines ; The program will first open the data file and read in data lines until the 'ZCZC' string is encountered. ; From that point it reads all lines in until the line with 'NNNN'. ; There are 24 'normal' comment lines generated from the header data, including ground check lines. NNCOML = 25 NCOML = strarr(nncoml) NCOML(0) = '+++++ Normal comments including MW11 header data ++++++' NCOML(1) = ' ' NCOML(2) = ' input filename: ' + txtfile NCOML(3) = ' # of data records: ' NCOML(4) = ' # of coded tropopauses: ' NCOML(5) = ' ' NCOML(6) = ' SOUNDING PROGRAM REV 8.35 USING GPS' NCOML(7) = '' NCOML(8) = ' STATION : 78762 ALAJUELA' NCOML(9) = ' LOCATION : 9.99 N 84.22 W 921 M' NCOML(10) = ' SOUNDING : ' ; with format = (i3) NCOML(11) = ' RS-NUMBER : ' ; with format = (a8) NCOML(12) = ' RADIOSONDE MODEL: ' ; with format = (a8) NCOML(13) = '' NCOML(14) = ' Ground Check : Ref RS Corr.' NCOML(15) = '' NCOML(16) = ' Pressure :' ; with format = (f11.1,f12.1,f9.1) NCOML(17) = ' Temperature :' ; with format = (f11.1,f12.1,f9.1) NCOML(18) = ' Humidity :' ; with format = (i11, i11, i9) NCOML(19) = '' NCOML(20) = ' STARTED AT : ' ; with ISO format YYYY-MM-DDTHHMM UTC NCOML(21) ='' NCOML(22) = ' Time AscRate Hgt/MSL Pressure Temp RH Dewp Dir Speed 88 CP 77 Theta U V Lat Lng QC flags' NCOML(23) = ' sec m/s m hPa degC % degC deg m/s K m/s m/s degN degE AHPTrdDS' NCOML(24) = '' XNAME = 'Seconds elapsed since sonde release (see "STARTED AT" below)' ; Primary variables NV = 24 VNAME = strarr(nv) VSCAL = fltarr(nv) VMISS = fltarr(nv) ; Basic variables VNAME(0) = 'Ascent rate (m/s)' VNAME(1) = 'Hgt/MSL (m)' VNAME(2) = 'Pressure (hPa)' VNAME(3) = 'Temperature (deg C)' VNAME(4) = 'Relative humidity (%)' VNAME(5) = 'Dewpoint (deg C)' VNAME(6) = 'Wind direction (deg)' VNAME(7) = 'Wind speed (m/s)' ; Profile characteristic points VNAME(8) = 'Coded tropopause (88 group) (0 = NO, 1 = first, 2 = second, 3 = 3rd)' VNAME(9) = 'Temp minimum (coldpoint) (0 = NO; 1 = YES)' VNAME(10) = 'Coded wind max (77 group) (0 = NO, 1 = YES)' ; Derived variables VNAME(11) = 'Potential temperature (K)' VNAME(12) = 'Zonal wind (derived) (m/s)' VNAME(13) = 'Meridional wind (derived) (m/s)' VNAME(14) = 'Latitude (derived) (deg N)' VNAME(15) = 'Longitude (derived) (deg E)' ; Quality control flags VNAME(16) = 'Quality flag - Ascent rate' VNAME(17) = 'Quality flag - Hgt/MSL' VNAME(18) = 'Quality flag - Pressure' VNAME(19) = 'Quality flag - Temperature' VNAME(20) = 'Quality flag - Relative humidity' VNAME(21) = 'Quality flag - Dewpoint' VNAME(22) = 'Quality flag - Wind direction' VNAME(23) = 'Quality flag - Wind speed' VSCAL = fltarr(nv) + 1.0 VMISS = fltarr(nv) VMISS = [ 99999.9, 999999, 99999.9, 9999.9, 999, 9999.9, 9999, 999.9, 99,99,99, $ 9999.99, 9999.99, 9999.99, 9999.999, 9999.999, 9, 9, 9, 9, 9, 9, 9, 9 ] ; -------------------------------------------------------------------------------------------------- ; ; READ THROUGH INPUT FILE ; ; -------------------------------------------------------------------------------------------------- ; First, read all the way through to the end of the file, which consists of three parts: ; (a) DigiCORA header data ; (b) 2-sec data ; (c) WMO message data ; In this part of the code, each of these parts is read into a separate array and processed accordingly inlines = strarr(4000) dataLines = strarr(4000) headLines = strarr(50) inline = '' iline = 0 nh = 0 N = 0 readHead = 1b readData = 0b readWMO = 0b readTTAA = 0b readTTCC = 0b readTTDD = 0b nTrop = 0 nWmax = 0 pressTrop = fltarr(4) tempTrop = fltarr(4) dTemp_prev = fltarr(4) + 10. NSCOML = 1 SCOML = strarr(100) SCOML(0) = '+++ Special comments including WMO FM-35 messages +++' SCOML(1) = ' ' zeroMinSec = ' 0 0' noRev = 1b irev = 0 comments = strarr(10) ncomments = 0 readLine: while not eof(12) do begin readf,12,inline iline = iline + 1 inlines(iline-1) = inline ; First, note the line in which is found the phrase "REV 8.3" is found ; Any non-blank lines ahead of this can be considered as comments to be added to the end of the special comment lines if noRev then begin lenline = strlen(inline) if lenline gt 0 then begin for j = 0, lenline-9 do begin if strupcase(strmid(inline,j,8)) eq 'REV 8.34' then begin NCOML(6) = ' SOUNDING PROGRAM REV 8.34 USING GPS' noRev = 0b iRev = iline if iRev gt 2 then goto, extractComments endif if strupcase(strmid(inline,j,8)) eq 'REV 8.35' then begin noRev = 0b iRev = iline if iRev gt 2 then goto, extractComments endif end endif extractComments: for i = 0, iRev-2 do begin ncomments = ncomments + 1 comments(i) = inlines(i) ; print,'Comment line ',i,':',inlines(i) end goto,readLine endif if strmid(inline,0,7) eq zeroMinSec then readData = 1b if strmid(inline,0,6) eq 'USCS01' then begin NSCOML = 2 readWMO = 1b readTTAA = 1b goto,readLine endif if readWMO then begin ; Check for tropopauses in TTAA messages if readTTAA then begin if strmid(inline,0,2) eq '88' and strmid(inline,2,3) ne '999' then begin nTrop = nTrop + 1 pressTrop(ntrop-1) = float(strmid(inline,2,3)) tempTrop(ntrop-1) = -float(strmid(inline,6,3)) / 10. endif if strmid(inline,12,2) eq '88' and strmid(inline,14,3) ne '999'then begin nTrop = nTrop + 1 pressTrop(ntrop-1) = float(strmid(inline,14,3)) tempTrop(ntrop-1) = -float(strmid(inline,18,3)) / 10. endif endif ; Check for tropopauses in TTCC messages if readTTCC then begin if strmid(inline,0,2) eq '88' and strmid(inline,2,3) ne '999' then begin nTrop = nTrop + 1 pressTrop(ntrop-1) = float(strmid(inline,2,3)) / 10. tempTrop(ntrop-1) = -float(strmid(inline,6,3)) / 10. endif if strmid(inline,12,2) eq '88' and strmid(inline,14,3) ne '999' then begin nTrop = nTrop + 1 pressTrop(ntrop-1) = float(strmid(inline,14,3)) / 10. tempTrop(ntrop-1) = -float(strmid(inline,18,3)) / 10. endif endif ; Check for windmax in TTAA messages if readTTAA then begin if strmid(inline,0,2) eq '77' and strmid(inline,2,3) ne '999' then begin nWmax = nWmax + 1 pressMax = float(strmid(inline,2,3)) wDirMax = fix(strmid(inline,6,2)) * 10 wDir01 = fix(strmid(inline,8,1)) if wDir01 eq 5 then wDirMax = wDirMax + 5 wSpdmax = fix(strmid(inline,8,3)) if wSpdMax gt 500 then wSpdMax = wSpdMax - 500 endif endif ; Check for windmax in TTCC messages if nWmax eq 0 and readTTCC then begin if strmid(inline,0,2) eq '77' and strmid(inline,2,3) ne '999' then begin nWmax = nWmax + 1 pressMax = float(strmid(inline,2,3)) / 10. wDirMax = fix(strmid(inline,6,2)) * 10 wDir01 = fix(strmid(inline,8,1)) if wDir01 eq 5 then wDirMax = wDirMax + 5 wSpdmax = fix(strmid(inline,8,3)) if wSpdMax gt 500 then wSpdMax = wSpdMax - 500 endif endif if readTTDD and strlen(inline) eq 0 then begin goto, readLine endif else if strmid(inline,0,6) eq 'UKCS01' then begin readTTAA = 0b goto, readLine endif else if strmid(inline,0,6) eq 'ULCS01' then begin readTTCC = 1b goto, readLine endif else if strmid(inline,0,6) eq 'UECS01' then begin readTTCC = 0b goto, readLine endif else if strmid(inline,6,4) eq 'TTDD' then begin readTTDD = 1b NSCOML = NSCOML + 1 SCOML(NSCOML-1) = inline goto,readLine endif else if strmid(inline,0,4) eq 'NNNN' then begin readWMO = 0b goto,readLine endif NSCOML = NSCOML + 1 SCOML(NSCOML-1) = inline ;print,SCOML(NSCOML-1) goto,readLine endif if readData then begin if strlen(inline) gt 0 then begin if strmid(inline,0,4) ne 'ZCZC' then begin dataLines(N) = inline ;if N le 10 then print,N,inline N = N + 1 endif else if N gt 1 then begin print,'ZCZC found. End of data after',N,' lines of data' ;for i = N-10, N-1 do print,i,datalines(i) endif endif goto, readLine endif if readHead then begin headLines(nh) = inline print,headLines(nh) nh = nh + 1 goto, readLine endif end ; Add prepended comments to the end of SCOML; otherwise just add a blank line if ncomments lt 2 then begin NSCOML = NSCOML + 1 SCOML(NSCOML-1) = '' endif else begin nextra = ncomments + 2 NSCOML = NSCOML + nextra SCOML(NSCOML-nextra-1) = '' SCOML(NSCOML-nextra) = '----------------- Notes/commentary ------------------' SCOML(NSCOML-nextra+1) = '' for j = 0, ncomments - 1 do SCOML(NSCOML-nextra+2+j) = comments(j) SCOML(NSCOML-1) = '' endelse ; Store the number of data lines in NCOML(3) ; Store the number of coded trops in NCOML(4) NCOML(3) = NCOML(3) + string(N) NCOML(4) = NCOML(4) + string(nTrop) close,12 ; -------------------------------------------------------------------------------------------------- ; ; DONE READING THROUGH INPUT FILE ; ; -------------------------------------------------------------------------------------------------- ; Results from tropopause and wind max decoding if nTrop gt 0 then begin print, 'Coded tropopauses: ' for i = 0, nTrop-1 do print,i,pressTrop(i),tempTrop(i) endif else print,' No coded tropopause (88999)' if nWmax gt 0 then begin print,' Wind max at ', pressMax,' dir/spd = ',wDirMax,wSpdMax endif else print,' No coded wind max (77999)' ; -------------------------------------------------------------------------------------------------- ; ; PROCESS DIGICORA HEADER DATA ; ; -------------------------------------------------------------------------------------------------- ; Extract sounding number, radiosonde serial and model numbers, and release date & time ih = 2 ; skip these two lines readSoundingNumber: while ih le nh do begin for j = 0, 5 do if strmid(strupcase(headLines(ih)),j,3) eq 'SOU' then goto, extractSoundingNumber ih = ih +1 goto, readSoundingNumber end print,'No sounding number found.' ih = 0 goto, readSerialNumber extractSoundingNumber: trimmedLine = strtrim(headLines(ih)) j0 = strlen(trimmedLine) - 4 soundingNumber = strtrim(strmid(trimmedLine,j0,4),1) NCOML(10) = NCOML(10) + soundingNumber print,NCOML(10) ih = 0 readSerialNumber: while ih le nh do begin for j = 0, 5 do $ if strmid(strupcase(headLines(ih)),j,3) eq 'RS-' then goto, extractSerialNumber ih = ih +1 goto, readSerialNumber end print,'No serial number found.' ih = 0 goto, readModelNumber extractSerialNumber: trimmedLine = strtrim(headLines(ih)) j0 = strlen(trimmedLine) - 9 serialNumber = strtrim(strmid(trimmedLine,j0,9),2) if strlen(serialNumber) ne 8 then begin print,'Extra character in radiosonde serial number: ',serialNumber endif NCOML(11) = NCOML(11) + serialNumber print,NCOML(11) ih = 0 readModelNumber: while ih le nh do begin for j = 0, 5 do $ if strmid(strupcase(headLines(ih)),j,3) eq 'RAD' then goto, extractModelNumber ih = ih +1 goto, readModelNumber end print,'No serial number found.' ih = 0 goto,readReleaseUTC extractModelNumber: trimmedLine = strtrim(headLines(ih)) j0 = strlen(trimmedLine) - 9 serialNumber = strtrim(strmid(trimmedLine,j0,9),2) NCOML(12) = NCOML(12) + serialNumber print,NCOML(12) ih = 0 readReleaseUTC: while ih le nh do begin for j = 0, 10 do $ if strmid(strupcase(headLines(ih)),j,5) eq 'START' then goto, extractReleaseUTC ih = ih + 1 goto, readReleaseUTC end print,'No release time found.' ih = 0 goto,readGroundCheck extractReleaseUTC: trimmedLine = strtrim(headLines(ih)) j0 = strlen(trimmedLine) - 25 releaseUTC = strtrim(strmid(trimmedLine,j0,25),2) print,'releaseUTC = ', releaseUTC ; Extract time from ReleaseUTC j0 = strlen(releaseUTC) - 9 releaseTime = strmid(releaseUTC,j0,5) print,'releaseTime = ',releaseTime hh = fix(strmid(releaseTime,0,2)) min = fix(strmid(releaseTime,3,2)) releaseHHMM = strcompress(string(hh) + string(min),/remove_all) if hh lt 10 then releaseHHMM = '0' + releaseHHMM dd = fix(strmid(releaseUTC,0,2)) releaseDD = strtrim(string(dd),2) if dd lt 10 then releaseDD = '0' + releaseDD print,'releaseDD = ',releaseDD releaseYYYY = '2005' if strmid(releaseUTC,2,1) eq ' ' then begin releaseMON = strupcase(strmid(releaseUTC,3,3)) endif else releaseMON = strupcase(strmid(releaseUTC,2,3)) ;print,'releaseMON = ', releaseMON, releaseUTC ;stop SMON = ['ENE','FEB','MAR','ABR','MAI','JUN','JUL','AGO','SET','OCT','NOV','DIC'] EMON = ['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEC'] mm = 0 oo = where( releaseMON eq SMON, nok ) if nok eq 1 then mm = oo(0) + 1 if nok ne 1 then oo = where( releaseMON eq EMON, nok) if nok eq 1 then mm = oo(0) + 1 if mm eq 0 then begin print, 'Bad releaseMON datum: ', releaseMON endif else begin releaseMM = string(mm) if mm lt 10 then releaseMM = '0' + releaseMM releaseUTC = releaseDD + ' ' + EMON(mm-1) + ' ' + releaseYYYY + ' ' + releaseTime + ' UTC' endelse NCOML(20) = NCOML(20) + releaseUTC ih = 0 readGroundCheck: ; The code used in the old version of 'edit_EDT.pro' to extract each of the three columns of ; data for each of the three ground check quantities, pressure,temperature, and humidity, ; was very long and complicated with lots of character manipulations. ; Here we simply look grab each of the lines with one of the unioque word fragments 'Press', ; 'Temp' and 'Humid' and deposit them into the appropriate comment lines. And be done with it! for ih = 2, nh-1 do begin for j = 0, 20 do begin if strmid(strupcase(headLines(ih)),j,5) eq 'PRESS' then begin hlen = strlen(headLines(ih)) NCOML(16) = strmid(headLines(ih),2,hlen-2) endif if strmid(strupcase(headLines(ih)),j,5) eq 'TEMPE' then begin hlen = strlen(headLines(ih)) NCOML(17) = ' T' + strmid(headLines(ih),6,hlen-6) endif if strmid(strupcase(headLines(ih)),j,5) eq 'HUMID' then begin hlen = strlen(headLines(ih)) NCOML(18) = strmid(headLines(ih),2,hlen-2) endif endfor endfor ; Open output file ---------------------------------- ; EDTfile = string(cdd) + string(cmm) + string(cyy) + string(chh) + 'z.EDT' ; EDTfile = strcompress(EDTfile,/remove_all) cyyyy = '20'+ cyy EDTfile = strcompress('RS78762_' + cyyyy + '-' + cmm + '-' + cdd + 'T' + chh + '.EDT',/remove_all) outPath = 'EDT_archive' + '/' + yyyyPathmm + '/' + hhZ + '/' + EDTfile print,outPath openw,22,outPath ; -------------------------------------------------------------------------------------------------- ; ; SET UP HEADER VARIABLES AND PRINT THEM TO OUTPUT FILE ; ; -------------------------------------------------------------------------------------------------- NHEADL = 9 + 2* fix((NV+0.5)/6) + NV + 1 + NSCOML + 1 + NNCOML ONAME = 'SELKIRK, HENRY' ORG = 'BAER INSTITUTE, M/S 245-5, NASA-AMES RES. CTR, MOFFETT FIELD, CA 94035-1000, (650) 604-6489, SNAME = 'VAISALA RADIOSONDE DATA - DigiCORA MW11 at Station MROC/78762, Alajuela, COSTA RICA; upgraded for RS92-SGP' MNAME = 'TICOSONDE-AURA/TCSP 2005' IVOL = 1 NVOL = 1 DATE = intarr(3) DATE(0) = fix(2000 + yy) DATE(1) = fix(mm) DATE(2) = fix(dd) RDATE = intarr(3) RDATE(0) = fix(strmid(systime(),20,4)) RMON = strupcase(strmid(systime(),4,3)) oo = where( EMON eq RMON, count) RDATE(1) = oo(0) + 1 RDATE(2) = fix(strmid(systime(),8,2)) DT = 2 printf,22,format='(2I6,T55,a)', NHEADL, 1001, '[ NHEADL = No. of header lines; FFI = Hipskind-Gaines format index ]' printf,22,format='(a,T55,a)', ONAME, '[ ONAME = Name of originator ]' printf,22,format='(a)', ORG printf,22,format='(a)', SNAME printf,22,format='(a,T55,a)', MNAME, '[ MISSION NAME ]' printf,22,format='(2I6,T55,a)', IVOL, NVOL, '[ IVOL NVOL ]' printf,22,format='(6I6,T55,a)', DATE,RDATE, '[ DATE, RDATE = Measurement date and processing date ]' printf,22,format='(f7.0,T55,a)', DT, '[ Dt (sec) = time base interval]' printf,22,format='(a,T55,a)', XNAME printf,22,format='(i6,T55,a)', NV, '[ NV = no. of original and derived variables/columns ]' printf,22, VSCAL printf,22, VMISS for i = 0, NV-1 do printf,22, VNAME(i) printf,22,format='(i6,T55,a)', NSCOML, '[ NSCOML = no. of special comment lines => WMO messages ]' for i = 0, NSCOML-1 do printf,22,SCOML(i) printf,22,format='(i6,T55,a)', NNCOML, '[ NNCOML = no. of normal comment lines => DigiCORA header data ]' for i = 0, NNCOML-1 do printf,22,NCOML(i) for i = 0, NSCOML-1 do print,SCOML(i) for i = 0, NNCOML-5 do print,NCOML(i) ; -------------------------------------------------------------------------------------------------- ; ; PROCESS SOUNDING DATA - N.B. STATION LAT and LON ARE HARDWIRED HERE ; ; -------------------------------------------------------------------------------------------------- seconds = 0 u0 = 0. v0 = 0. lat0 = 9.99 lon0 = -84.22 dXkm = 0. dYkm = 0. tempColdPt = 9999. HgtTOP = 0L mins0 = 0 secs0 = 0 ; Initialize the sounding variables to missing values TimeSec = lonarr(N) AscRate = fltarr(N) + VMISS(0) HgtMSL = lonarr(N) + VMISS(1) Press = fltarr(N) + VMISS(2) Temp = fltarr(N) + VMISS(3) RelHum = intarr(N) + VMISS(4) DewPt = fltarr(N) + VMISS(5) Wdir = intarr(N) + VMISS(6) Wspd = fltarr(N) + VMISS(7) ; Initialize the quality flags qAscRate = bytarr(N) + 0 qHgtMSL = bytarr(N) + 0 qPress = bytarr(N) + 0 qTemp = bytarr(N) + 0 qRelHum = bytarr(N) + 0 qDewPt = bytarr(N) + 0 qWdir = bytarr(N) + 0 qWspd = bytarr(N) + 0 gotWind = 0b i0 = -1 Theta = fltarr(N) + VMISS(11) U = fltarr(N) + VMISS(12) V = fltarr(N) + VMISS(13) Lat = fltarr(N) + VMISS(14) Lon = fltarr(N) + VMISS(15) kappa = 2./7. ; -------------------------------------------------------------------------------------------------- ; ; LOOP THROUGH 2-sec DATA, CHECK FOR MISSING DATA ; ; -------------------------------------------------------------------------------------------------- for i = 0, N-1 do begin mins = fix(strtrim(strmid(dataLines(i),0,4),2)) secs = fix(strtrim(strmid(dataLines(i),5,2),2)) strAscRate = strmid(dataLines(i), 7,8) strHgtMSL = strmid(dataLines(i),15,8) strPress = strmid(dataLines(i),23,9) strTemp = strmid(dataLines(i),32,7) strRelHum = strmid(dataLines(i),39,4) strDewPt = strmid(dataLines(i),43,7) strWdir = strmid(dataLines(i),50,5) strWspd = strmid(dataLines(i),55,6) ; Check for slash (/) characters which indicate missing data slash = 0b for j = 0,7 do if strmid(strAscRate,j,1) eq '/' then slash = 1b if not slash then AscRate(i) = float(strAscRate) slash = 0b for j = 0,7 do if strmid(strHgtMSL,j,1) eq '/' then slash = 1b if not slash then HgtMSL(i) = long(strHgtMSL) slash = 0b for j = 0,8 do if strmid(strPress,j,1) eq '/' then slash = 1b if not slash then Press(i) = float(strPress) slash = 0b for j = 0,6 do if strmid(strTemp,j,1) eq '/' then slash = 1b if not slash then Temp(i) = float(strTemp) slash = 0b for j = 0,3 do if strmid(strRelHum,j,1) eq '/' then slash = 1b if not slash then RelHum(i) = fix(strRelHum) slash = 0b for j = 0,6 do if strmid(strDewPt,j,1) eq '/' then slash = 1b if not slash then DewPt(i) = float(strDewPt) slash = 0b for j = 0,4 do if strmid(strWdir,j,1) eq '/' then slash = 1b if not slash then Wdir(i) = fix(strWdir) slash = 0b for j = 0,5 do if strmid(strWspd,j,1) eq '/' then slash = 1b if not slash then Wspd(i) = float(strWspd) ; -------------------------------------------------------------------------------------------------- ; ; DATA QUALITY CHECKS ; ; -------------------------------------------------------------------------------------------------- if AscRate(i) lt 0. then begin qAscRate(i) = 1b print, mins, secs, AscRate(i),' - AscRate < 0 - Set quality flag' endif if HgtMSL(i) gt 50000 and HgtMSL(i) ne VMISS(1) then begin qHgtMSL(i) = 1b print, mins, secs, HgtMSL(i),' - HgtMSL > 50 km - Set quality flag and reset HgtMSL to missing value' HgtMSL(i) = VMISS(1) endif else if HgtMSL(i) lt 0 then begin qHgtMSL(i) = 1b print, mins, secs, HgtMSL(i),' - HgtMSL < 0 - Set quality flag and reset HgtMSL to missing value' HgtMSL(i) = VMISS(1) endif if Press(i) gt 1000. and Press(i) ne VMISS(2) then begin qPress(i) = 1b print, mins, secs, Press(i),' - Press > 1000 hPa - Set quality flag and reset Press to missing value' Press(i) = VMISS(2) endif else if Press(i) lt 0. then begin qPress(i) = 1b print, mins, secs, Press(i) ,' - Press < 0 - Set quality flag and reset Press to missing value' Press(i) = VMISS(2) endif if Temp(i) gt 50. and Temp(i) ne VMISS(3) then begin qTemp(i) = 1b print, mins, secs, Temp(i) ,' - Temp > 50 degC - Set quality flag and set Temp to missing value' Temp(i) = VMISS(3) endif else if Press(i) gt 800. and Temp(i) le 10. then begin qTemp(i) = 1b print, mins, secs, Temp(i) ,' - Press > 800 and Temp < 10 degC - Set quality flag' endif else if Press(i) gt 600. and Temp(i) le -10. then begin qTemp(i) = 1b print, mins, secs, Temp(i) ,' - Press > 600 and Temp < -10 degC - Set quality flag' endif else if Press(i) gt 400. and Temp(i) le -25. then begin qTemp(i) = 1b print, mins, secs, Temp(i) ,' - Press > 400 and Temp < -25 degC - Set quality flag' endif else if Press(i) gt 200. and Temp(i) le -65. then begin qTemp(i) = 1b print, mins, secs, Tem(i) ,' - Press > 200 and Temp < -65 degC - Set quality flag' endif else if Temp(i) lt -100. then begin qTemp (i) = 1b print, mins, secs, Temp(i) ,' - Temp < -100 degC - Set quality flag' endif if RelHum(i) lt 0 then begin qRelHum(i) = 1b print, mins, secs, RelHum(i) ,' - RelHum < 0 - Set quality flag and reset RelHum to missing value' RelHum(i) = VMISS(4) endif if Dewpt(i) gt 50. and DewPt(i) ne VMISS(5) then begin qDewPt(i) = 1b print, mins, secs, Dewpt(i) ,' - Dewpt > 50 degC - Set quality flag and reset DewPt to missing value' DewPt(i) = VMISS(5) endif if Wdir(i) gt 360 and Wdir(i) ne VMISS(6) then begin qWdir(i) = 1b print, mins, secs, Wdir(i) ,' - Wdir > 360 - Set quality flag and reset Wdir to missing value' Wdir(i) = VMISS(6) endif else if WDir(i) lt 0 then begin qWdir(i) = 1b print, itime,WDir(i) ,' - WDir < 0 - Set quality flag and reset Wdir to missing value' Wdir(i) = VMISS(6) endif if WSpd(i) gt 100. and Wspd(i) ne VMISS(7) then begin qWspd(i) = 1b print, seconds, WSpd(i) ,' - WSpd > 100 m/s - Set quality flag' endif else if WSpd(i) lt 0. then begin qWspd(i) = 1b print, seconds, WSpd(i) ,' - WSpd < 0 - Set quality flag and reset Wspd to missing value' Wspd(i) = VMISS(7) endif ; -------------------------------------------------------------------------------------------------- ; ; COMPUTE POTENTIAL TEMPERATURE ; ; -------------------------------------------------------------------------------------------------- if Press(i) ne VMISS(2) and Temp(i) ne VMISS(3) then begin Theta(i) = (Temp(i) +273.15) * (1000./Press(i))^kappa endif else Theta(i) = VMISS(11) ; -------------------------------------------------------------------------------------------------- ; ; INTEGRATE LAT AND LON IF WE HAVE WIND SOMEWHERE BELOW 900 mb ; ; -------------------------------------------------------------------------------------------------- ; if Wdir(i) eq VMISS(6) or Wspd(i) eq VMISS(7) or qWdir(i) or qWspd(i) then begin ; print, i, Wdir(i), Wspd(i), qWdir(i), qWspd(i) ; ;stop ; goto, nextLevel ; endif TimeSec(i) = 60L * mins + secs if Wdir(i) eq VMISS(6) or Wspd(i) eq VMISS(7) or qWdir(i) or qWspd(i) then begin ;print, TimeSec(i), Wdir(i), Wspd(i), qWdir(i), qWspd(i) ;stop goto, nextLevel endif if not gotWind and mins lt 5 then begin print,'Found first wind level: ',i,mins,secs,TimeSec(i),HgtMSL(i),Press(i),Wdir(i),Wspd(i) gotWind = 1b ;stop endif if i gt 0 and i0 ge 0 and gotWind then begin dSec = TimeSec(i) - TimeSec(i0) if dSec gt Dt then begin print,'Wind data gap begins at t0 =',TimeSec(i0+1),' Z0=',fix(HgtMSL(i0+1)),' dSec =',dSec endif else if dSec lt 0 then begin print,'Wind data repeat: i0=',i0,' t0=',TimeSec(i0),' Z0=',HgtMSL(i0) $ ,'i1=',i ,' t1=',TimeSec(i), ' Z1=',HgtMSL(i) endif endif else dSec = 0 mins0 = mins secs0 = secs angle = Wdir(i) + 90. rangle = angle*3.14159265358979/180. u(i) = Wspd(i) * cos(rangle) v(i) = -Wspd(i) * sin(rangle) maxGap = 600 if gotWind and dSec gt maxGap then begin print,'Gap =',dSec,' > ',maxGap,' so terminate lat/lon integration.' gotWind = 0b endif if not gotWind then goto, nextLevel i0 = i ; The index to the previous level with wind ; This set AFTER we read in the first good wind level and ; find u and v at that level. At the first good wind level ; dSec = 0, so no integration occurs until the second goo ; good wind level is found. if TimeSec(i) ne 0 then begin dXkm = dXkm + (u(i0)+u(i))/2 * dSec / 1000. dYkm = dYkm + (v(i0)+v(i))/2 * dSec / 1000. lat(i) = lat0 + dYkm / 1.11e2 rlat = lat(i) * 3.14159265358979 / 180. lon(i) = lon0 + dXkm / (1.11e2 * cos(rlat)) endif else begin lat(i) = lat0 lon(i) = lon0 endelse nextLevel: end ; -------------------------------------------------------------------------------------------------- ; ; SECOND TIME THROUGH THE DATA ; ; LOOK FOR CODED TROP AND WMAX IN 2-sec DATA; ALSO SEARCH FOR TEMP MINIMUM ; ; -------------------------------------------------------------------------------------------------- ; Initilize the special level data trop88 = intarr(N) + 0 coldPt = intarr(N) + 0 wMax77 = intarr(N) + 0 iTrop = 0 ilast = 0 k = 1b imax = -1 for i = 0, N-1 do begin checkTrop: if nTrop gt 0 and iTrop le nTrop then begin dPress = abs(Press(i) - pressTrop(itrop)) if Press(i) ge 100. and dPress lt 1 then begin if k then print,'Searching for trop # ',itrop+1,' coded as ',pressTrop(itrop),tempTrop(itrop) k = 0b dTemp = abs( Temp(i) - tempTrop(itrop) ) ;print,'dTemp = ', dTemp ; Case 1 - Temp(i) is odd - definite trop if abs(dTemp) lt 0.05 then begin print,'Definite trop ', itrop + 1,' at ',Press(i),' hPa and ',Temp(i),' deg C' tempTrop(itrop) = Temp(i) pressTrop(itrop) = Press(i) trop88(i) = iTrop + 1 iTrop = iTrop + 1 k = 1b ; Case 2 - Temp(i) is even and 0.1 deg warmer than tempTrop - need to check for next level too endif else begin print,'Possible trop ', itrop + 1,' at i=',i,Press(i),' hPa and ',Temp(i),' deg C.' dTempMin = dTemp iiTempMin = i for ii = i+1, i+10 do begin ; look for any of the next ten levels above for closer Temp dTemp1 = abs(Temp(ii) - tempTrop(itrop)) ;print,ii,Press(ii),Temp(ii),dTemp1,dTemp if dTemp1 lt dTemp and abs(Press(ii)-pressTrop(itrop)) lt 1 then goto, checkWmax end print,'OK, nothing better up to ',Press(ii),' so trop ',itrop+1,' at ',Press(i),' is definite.' tempTrop(itrop) = Temp(i) pressTrop(itrop) = Press(i) trop88(i) = iTrop + 1 iTrop = iTrop + 1 k = 1b endelse endif else if dPress eq 0.0 then begin ; if k then print,'Searching for trop # ',itrop+1,' coded as ',pressTrop(itrop),tempTrop(itrop) print,'Definite trop ', itrop + 1,' at ',Press(i),' hPa and ',Temp(i),' deg C tempTrop(itrop) = Temp(i) pressTrop(itrop) = Press(i) trop88(i) = iTrop + 1 iTrop = iTrop + 1 k = 1b endif endif checkWmax: if nWmax eq 0 then goto, checkColdPt dPress = abs(Press(i) - pressMax) if Press(i) ge 100 then begin if dPress lt 1 then begin if ilast gt 0 then begin if Wspd(i) gt Wspd(ilast) and Wspd(i) ne VMISS(7) then begin wMax77(ilast) = 0 wMax77(i) = 1 ilast = i imax = i print,'Coded wMax: ',pressMax, wDirMax, wSpdMax, 'kts' print,'Next wMax in 2-sec data: ', Press(i), Wdir(i), Wspd(i)*1.93, ' kts' endif endif else begin wMax77(i) = 1 ilast = i imax = i print,'Coded wMax: ',pressMax, wDirMax, wSpdMax, 'kts' print,'First wMax in 2-sec data: ', Press(i), Wdir(i), Wspd(i)*1.93, ' kts' endelse endif endif else begin if dPress eq 0 then begin if ilast gt 0 then begin if Wspd(i) gt Wspd(ilast) and Wspd(i) ne VMISS(7) then begin wMax77(ilast) = 0 wMax77(i) = 1 ilast = i imax = i print,'Coded wMax: ',pressMax, wDirMax, wSpdMax, 'kts' print,'Next wMax in 2-sec data: ', Press(i), Wdir(i), Wspd(i)*1.93, ' kts' endif endif else begin wMax77(i) = 1 ilast = i imax = i print,'Coded wMax: ',pressMax, wDirMax, wSpdMax, 'kts' print,'First wMax in 2-sec data: ', Press(i), Wdir(i), Wspd(i)*1.93, ' kts' endelse endif endelse checkColdPt: ; Check if this is the min temp so far if Temp(i) lt tempColdPt then begin tempColdPt = Temp(i) pressColdPt = Press(i) endif if HgtMSL(i) ne VMISS(1) then HgtTop = HgtMSL(i) end ; -------------------------------------------------------------------------------------------------- ; ; OUTPUT DATA ; ; -------------------------------------------------------------------------------------------------- if nwMax gt 0 and imax lt 0 then print,'Probable truncation of data below coded wind max of ' $ ,pressMax,'hPa ',wDirMax,' deg ', wSpdMax,' kts' for i = 0, N-1 do begin if Press(i) ne pressColdPt or pressColdPt ge 200. then begin coldPt(i) = 0 endif else begin if nTrop gt 0 then begin for j = 0, nTrop-1 do print,'Trop ',j+1,' at ', pressTrop(j), tempTrop(j) endif if HgtMSL(N-1)-HgtMSL(i) gt 1000 then begin print,'Cold point at ', Press(i), HgtMSL(i), Temp(i) coldPt(i) = 1 endif endelse if nwMax gt 0 then begin if i eq imax then print,'Wind max at ', imax, Press(i), wDir(i), wSpd(i), ' m/s' endif if HgtMSL(i) eq HgtTOP then print,'Maximum altitude in 2-sec data: ', HgtTOP, ' at P = ', Press(i) printf,22, format='(i7,f8.1,i8,f9.1,f7.1,i4,f7.1,i5,f6.1,2x,3i3,2x,3f9.2,2x,2f10.3,2x,8i1 )' $ , TimeSec(i), AscRate(i), HgtMSL(i), Press(i), Temp(i), RelHum(i), DewPt(i), Wdir(i), Wspd(i) $ , Trop88(i), ColdPt(i), Wmax77(i), Theta(i), U(i), V(i), Lat(i), Lon(i) $ , qAscRate(i), qHgtMSL(i), qPress(i), qTemp(i), qRelHum(i), qDewPt(i), qWdir(i), qWspd(i) end close, 22 end