/********************************************************************************* dstlsonic.mac FAN script file for viewing/processing DSTL sonic data files Kosteniuk Consulting Ltd Jan., 2002 Remarks: The DSTL sonic data used by this script is stored in binary files (one file per day) with 32 channels (rows) and variable length (columns). The 20-Hz time series for each 3-component sonic anemometer occupies 4-channels: U,V,W,T The order of the (8) sonics in the binary data files is: A1 (2 meter) A2 (6 meter) B1 (2 meter) B2 (6 meter) C1 (2 meter) C2 (6 meter) D1 (2 meter) D2 (6 meter) where Ai refers to the northwest tower Bi refers to the southwest tower Ci refers to the northeast tower Di refers to the southeast tower *********************************************************************************/ /*********************************************************** load the graphics/printer layout file (landscape version) ***********************************************************/ run {"graphicsl.mac"} /*********************************************************** define the directories containing the MUST data ***********************************************************/ MROOT = "/home/dpid/must1/" /* root directory for all MUST data */ DSTSD = MROOT + {"dstlsonic/"} /* directory containing DSTL sonic files*/ DSTBD = MROOT + {"dstlsonic/bin/"} /* directory containing binary files */ DSTOD = MROOT + {"dstlsonic/stats/"} /* directory containing statistic files */ DSTPD = MROOT + {"dstlsonic/plots/"} /* directory containing plot files */ DSTAD = MROOT + {"dstlsonic/ascii/"} /* directory containing ascii data files*/ /*********************************************************** define the names of the binary DSTL data files (1 file/day) ***********************************************************/ SFILES = "m2562043\n" /* Sept 13 */ SFILES = SFILES + {"m2580615\n"} /* Sept 15 */ SFILES = SFILES + {"m2610136\n"} /* Sept 18 */ SFILES = SFILES + {"m2620130\n"} /* Sept 19 */ SFILES = SFILES + {"m2630226\n"} /* Sept 20 */ SFILES = SFILES + {"m2640116\n"} /* Sept 21 */ SFILES = SFILES + {"m2650215\n"} /* Sept 22 */ SFILES = SFILES + {"m2671840\n"} /* Sept 24 */ SFILES = SFILES + {"m2681821\n"} /* Sept 25 */ SFILES = SFILES + {"m2692015"} /* Sept 26 */ /*********************************************************** define the start date/time of the binary DPD files ***********************************************************/ SDATES = [2001 256 204300 ] \ [2001 258 061500 ] \ [2001 261 013600 ] \ [2001 262 013000 ] \ [2001 263 022600 ] \ [2001 264 011600 ] \ [2001 265 021500 ] \ [2001 267 184000 ] \ [2001 268 182100 ] \ [2001 269 201500 ] /*********************************************************** define the rotation angles (in degrees cw from north) and height above ground level of each 2- or 3-component sonic ***********************************************************/ SRAHTS = [178.1 2.0 ] /* A1 */ \ [177.0 6.0 ] /* A2 */ \ [184.7 2.0 ] /* B1 */ \ [184.0 6.0 ] /* B2 */ \ [179.5 2.0 ] /* C1 */ \ [178.0 6.0 ] /* C2 */ \ [186.3 2.0 ] /* D1 */ \ [185.0 6.0 ] /* D2 */ /*********************************************************** define the names of (used in plots ...) of each sonic ***********************************************************/ SNAMES = "A1\n" /* NW 2 m */ SNAMES = SNAMES + {"A2\n"} /* NW 6 m */ SNAMES = SNAMES + {"B1\n"} /* SW 2 m */ SNAMES = SNAMES + {"B2\n"} /* SW 6 m */ SNAMES = SNAMES + {"C1\n"} /* NE 2 m */ SNAMES = SNAMES + {"C2\n"} /* NE 6 m */ SNAMES = SNAMES + {"D1\n"} /* SE 2 m */ SNAMES = SNAMES + {"D2\n"} /* SE 6 m */ /*********************************************************** define the date/times/... of the sonic data to process col 1 start year (yyyy) of the sonic data file col 2 start day (ddd) of the sonic data file col 3 start time (hhmmss) of the sonic data file col 4 start time (hhmmss) of the dpid/uvic data file col 5 offset (seconds relative to col 4) to start processing col 6 offset (seconds relative to col 4) to end processing col 7 postfix (1:?) of the output sonic stats file ***********************************************************/ PDATES = [2001 256 204300 221000 200 600 1 ] \ [2001 256 204300 221000 700 1100 2 ] \ [2001 256 204300 221000 300 500 3 ] \ [2001 256 204300 221000 800 1000 4 ] \ [2001 256 204300 223500 200 1000 1 ] \ [2001 256 204300 223500 300 500 2 ] \ [2001 256 204300 223500 500 700 3 ] \ [2001 256 204300 225700 200 700 1 ] \ [2001 256 204300 225700 700 1000 2 ] \ [2001 256 204300 225700 300 500 3 ] \ [2001 256 204300 232500 500 1000 1 ] \ [2001 256 204300 232500 700 900 2 ] \ [2001 261 013600 054300 1000 1200 1 ] \ [2001 261 013600 073100 300 500 1 ] \ [2001 261 013600 073100 800 1000 2 ] \ [2001 261 013600 075800 200 400 1 ] \ [2001 262 013000 022400 200 650 1 ] \ [2001 262 013000 022400 800 1000 2 ] \ [2001 262 013000 022400 300 500 3 ] \ [2001 262 013000 024600 200 700 1 ] \ [2001 262 013000 024600 500 700 2 ] \ [2001 262 013000 030700 500 900 1 ] \ [2001 262 013000 030700 600 800 2 ] \ [2001 262 013000 035600 600 900 1 ] \ [2001 262 013000 035600 650 850 2 ] \ [2001 263 022600 070100 800 1200 1 ] \ [2001 263 022600 070100 800 1000 1 ] \ [2001 264 011600 013800 500 700 1 ] \ [2001 264 011600 024600 450 800 1 ] \ [2001 264 011600 024600 1100 1300 2 ] \ [2001 264 011600 024600 500 700 3 ] \ [2001 264 011600 070100 300 500 1 ] \ [2001 264 011600 070100 400 600 2 ] \ [2001 267 184000 185200 200 1300 1 ] \ [2001 267 184000 185200 400 600 2 ] \ [2001 267 184000 193400 400 900 1 ] \ [2001 267 184000 193400 400 600 2 ] \ [2001 267 184000 200300 400 600 1 ] \ [2001 267 184000 203300 400 900 1 ] \ [2001 267 184000 203300 400 600 2 ] \ [2001 267 184000 210100 300 600 1 ] \ [2001 267 184000 210100 350 550 2 ] \ [2001 267 184000 215000 150 1000 1 ] \ [2001 267 184000 215000 150 500 2 ] \ [2001 267 184000 215000 300 500 3 ] \ [2001 267 184000 221300 75 350 1 ] \ [2001 267 184000 221300 600 900 2 ] \ [2001 267 184000 221300 100 300 3 ] \ [2001 267 184000 223500 50 400 1 ] \ [2001 267 184000 223500 400 600 2 ] \ [2001 267 184000 223500 200 400 3 ] \ [2001 267 184000 230300 300 500 1 ] \ [2001 267 184000 230300 700 900 2 ] \ [2001 268 182100 182900 100 900 1 ] \ [2001 268 182100 182900 300 500 2 ] \ [2001 268 182100 184900 50 950 1 ] \ [2001 268 182100 184900 100 300 2 ] \ [2001 268 182100 215000 150 950 1 ] \ [2001 268 182100 215000 300 500 2 ] \ [2001 268 182100 221100 200 900 1 ] \ [2001 268 182100 221100 500 700 2 ] \ [2001 268 182100 225600 50 900 1 ] \ [2001 268 182100 225600 300 500 2 ] \ [2001 268 182100 232000 150 950 1 ] \ [2001 268 182100 232000 500 700 2 ] \ [2001 268 182100 235300 150 950 1 ] \ [2001 268 182100 235300 400 600 2 ] \ [2001 269 201500 205400 400 900 1 ] \ [2001 269 201500 205400 1000 1200 2 ] \ [2001 269 201500 213100 150 1100 1 ] \ [2001 269 201500 213100 500 700 2 ] \ [2001 269 201500 215700 200 900 1 ] \ [2001 269 201500 215700 400 600 2 ] \ [2001 269 201500 222300 200 400 1 ] \ [2001 269 201500 222300 700 900 2 ] \ [2001 269 201500 225000 200 800 1 ] \ [2001 269 201500 225000 300 500 2 ] \ /*********************************************************** define some constants used in the sonic data processing ***********************************************************/ DSRATE = 20 /* DSTL sonic sampling rate (Hz) */ NSONIC = 8 /* total number of DSTL sonics */ GNORTH = 30.0 ; /* angle (deg) between grid and true north */ GNORTH = 0.0 ; /* needed - rotation angles are relative to grid*/ PI = 3.141592654 /* Pi */ GRAVC = 9.81 /* gravitational constant (m/s^2) */ DTORAD = PI / {180.0} /* conversion factor for degrees to radians */ RADTOD = {180.0} / PI /* conversion factor for radians to degrees */ MWDA = 28.9666 /* molecular weight of dry air (kg/kmol) */ GAMMA = 1.40 /* ratio of specific heat capacities of dry air */ FCONST = 0.9997 /* dimensionless constant */ VONKC = 0.4 /* von Karman's constant */ DENAIR = 1.1 /* density of air at DPG (kg/m^3) */ SHCAIR = 1005.0 /* specific heat cap of air (J/kg-K) */ RCONST = 8314.36 /* universal gas constant (J/kmol-K) */ /* constant to convert speed of sound to temp */ SOSTOK = MWDA / {GAMMA * {FCONST * RCONST}} NBLOCK = 600 /* block average to 1 sample / 30 s for plots */ /*********************************************************** Part 2: Define the macros (functions) accessed from the menu ***********************************************************/ /******************************************************************************* loadbfile - Load a binary sonic data file Inputs: SCODE index (0:?) of the sonic file to load Outputs: sdata 52-row matrix containing the binary sonic data sdate the start time (in unix calander format) *******************************************************************************/ macro loadbfile ; \ SFILE = extract SFILES SCODE ; /* base name of the data file */ \ SDATE = extract SDATES SCODE ; /* start date/time (y d hhmmss) */ \ sfile = DSTBD + {SFILE + {".bin"}}; /* full name of the data file */ \ rfile = DSTBD + {SFILE + {".cor"}}; /* full name for corrected data */ \ echo {" "} ; /* output a blank line */ \ echo {{"reading file "}+sfile} ; /* send the user a message */ \ sdata = bread sfile ; /* read the binary file */ \ sdate = mktime SDATE ; /* convert date to cal format */ \ panel pid ; /* redraw the panel */ /******************************************************************************* rotabfile - Rotate the sonic data file *******************************************************************************/ macro rotabfile ; \ rdata = sdata ; /* make a copy of the data */ \ nsamp = col sdata ; /* number of samples/channel */ \ nsonic = NSONIC ; while nsonic < /* loop over sonics */ \ isonic = NSONIC - nsonic ; /* sonic index (0:NSONIC-1) */ \ sanght = extract SRAHTS isonic ; /* get the sonic angle/height */ \ sangle = extract sanght {0} ; /* get the rotation angle (deg) */ \ sheight= extract sanght {1} ; /* get the sonic height (m AGL) */ \ sangle = sangle + GNORTH ; /* rotation angle wrt grid north*/ \ rangle = sangle - {180.0} ; /* angle from +U axis to GN */ \ rangle = rangle * DTORAD ; /* convert from degrees to rads */ \ sname = extract SNAMES isonic ; /* get the name of the sonic */ \ r1 = {4}*isonic ; /* first row to extract */ \ r2 = r1 + {3} ; /* last row to extract */ \ c1 = 0 ; /* first col to extract */ \ c2 = nsamp - {1} ; /* last col to extract */ \ rc1 = vector r1 c1 ; /* combine first row/col */ \ rc2 = vector r2 c2 ; /* combine last row/col */ \ tdata = extract sdata rc1 rc2 ; /* extract current sonic data */ \ ucomp = {.01} * {extract tdata {0}}; /* u component (m/s) */ \ vcomp = {.01} * {extract tdata {1}}; /* v component (m/s) */ \ wcomp = {.01} * {extract tdata {2}}; /* w component (m/s) */ \ ccomp = {.01} * {extract tdata {3}}; /* c component (m/s) */ \ sra = sin rangle ; /* sin(rotation angle) (-1:1) */ \ cra = cos rangle ; /* cos(rotation angle) (-1:1) */ \ gns = {cra*ucomp} - {sra*vcomp} ; /* wind speed along Gnorth axis */ \ ges = {sra*ucomp} + {cra*vcomp} ; /* wind speed along Geast axis */ \ gws = {-1.0} * ges ; /* wind speed along Gwest axis */ \ stm = SOSTOK * {ccomp*ccomp} ; /* convert sos to sonic temp (K)*/ \ replace gns tdata {0 0} ; /* save wind speed along Gnorth */ \ replace gws tdata {1 0} ; /* save wind speed along Gwest */ \ replace wcomp tdata {2 0} ; /* save vertical wind speed */ \ replace stm tdata {3 0} ; /* save wind speed along Gwest */ \ replace tdata rdata rc1 ; /* save correct sonic data */ \ nsonic = nsonic - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next sonic */ \ echo {" "} ; /* output a blank line */ \ echo {{"creating file "}+rfile} ; /* send the user a message */ \ bwrite rdata rfile ; /* write the corrected data file*/ \ panel pid ; /* redraw the panel */ /************************************************************************* compbstat - Compute the basic stats for a single trial/period Inputs: PCODE index (0:?) of the trial/period to process *************************************************************************/ macro compbstat ; \ PDATE = extract PDATES PCODE ; /* get the trial date/time info */ \ syyy = extract PDATE {0} ; /* sonic file year */ \ sddd = extract PDATE {1} ; /* sonic file day */ \ shms = extract PDATE {2} ; /* sonic file hour/min/sec */ \ thms = extract PDATE {3} ; /* trial hour/min/sec */ \ toff = extract PDATE {4} ; /* start offset to process */ \ tend = extract PDATE {5} ; /* end offset to process */ \ tnum = extract PDATE {6} ; /* postfix of the output file */ \ tdur = tend - toff ; /* length (s) of the series */ \ /* convert to a sonic filename */ \ sdhms = {sddd*{10000}}+{shms/{100}}; /* dddhhmm */ \ rfile = {{"m"}+{ascii sdhms}}+{".cor"}; /* sonic file = mdddhhmm.cor */ \ /* convert to an stats filename */ \ odhms = {sddd*{10000}}+{thms/{100}}; /* dddhhmm */ \ ofile = {{"m"}+{ascii odhms}}+{"_"}; /* stats file = mdddhhmm */ \ ofile = ofile+{{ascii tnum}+{".txt"}};/* stats file = mdddhhmm_n.txt */ \ rfile = DSTBD + rfile ; /* full name of the sonic file */ \ ofile = DSTOD + ofile ; /* full name of the stats file */ \ /* determine the starting offset*/ \ scal = mktime{vector syyy sddd shms};/* start of sonic (cal format) */ \ tcal = mktime{vector syyy sddd thms};/* start of trial (cal format) */ \ soff = {tcal - scal} + toff ; /* offset (s) from sonic start */ \ c1 = DSRATE * soff ; /* first column to extract */ \ c2 = c1 + {{DSRATE*tdur} - {1}} ; /* last column to extract */ \ r1 = 0 ; /* first row to extract */ \ r2 = {{4}*NSONIC} - {1} ; /* last row to extract */ \ rc1 = vector r1 c1 ; /* combine first row/col */ \ rc2 = vector r2 c2 ; /* combine last row/col */ \ echo {" "} ; /* output a blank line */ \ echo {{"loading sonic file "}+rfile};/* print the sonic file name */ \ rdata = bread rfile ; /* load the entire sonic file */ \ rdata = extract rdata rc1 rc2 ; /* extract the desired portion */ \ echo {{" processing "} + {{ascii {col rdata}} + {" samples"}}} ; \ echo {{"creating stats file "} + ofile};/* print the stats file name */ \ report ofile; /* create the stats file */ \ /* write the file header */ \ tname = {ascii odhms} + {"M"} ; /* trial name (like 2562210M) */ \ tperd = {ascii toff}+{{":"}+{ascii tend}} ; /* period */ \ print {{"Trial "}+tname} ; /* write the trial name */ \ print {{"Period "}+{tperd+{"s"}}} ; /* write the period */ \ print {" "} ; /* blank line */ \ /* assemble the column header */ \ thead = "Sonic " ; /* start with 8 spaces */ \ nsonic = NSONIC ; while nsonic < /* loop over sonics */ \ isonic = NSONIC - nsonic ; /* sonic index (0:NSONIC-1) */ \ soname = extract SNAMES isonic ; /* name of the sonic (N1,...) */ \ thead = thead + {{" "} + soname}; /* add the sonic name */ \ nsonic = nsonic - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next sonic */ \ print thead ; /* output the column header */ \ print {" "} ; /* blank line */ \ /* compute stats for all sonics */ \ tstat = matrix {44} NSONIC ; /* storage for results */ \ nsamp = col rdata ; /* number of samples/channel */ \ nsonic = NSONIC ; while nsonic < /* loop over sonics */ \ isonic = NSONIC - nsonic ; /* sonic index (0:NSONIC-1) */ \ sanght = extract SRAHTS isonic ; /* get the sonic angle/height */ \ sheight= extract sanght {1} ; /* get the sonic height (m AGL) */ \ r1 = {4}*isonic ; /* first row to extract */ \ r2 = r1 + {3} ; /* last row to extract */ \ c1 = 0 ; /* first col to extract */ \ c2 = nsamp - {1} ; /* last col to extract */ \ rc1 = vector r1 c1 ; /* combine first row/col */ \ rc2 = vector r2 c2 ; /* combine last row/col */ \ tdata = extract rdata rc1 rc2 ; /* extract current sonic data */ \ U = extract tdata {0}; /* x component (m/s) */ \ V = extract tdata {1}; /* y component (m/s) */ \ W = extract tdata {2}; /* z component (m/s) */ \ T = extract tdata {3}; /* temperature (K) */ \ A = RADTOD * {atan U V} ; /* bearing (deg CCW from Grid N)*/ \ /* compute means of U ... */ \ Ubar = mean U ; /* mean x component */ \ Vbar = mean V ; /* mean y component */ \ Wbar = mean W ; /* mean z component */ \ Tbar = mean T ; /* mean temperature */ \ Abar = RADTOD * {atan Ubar Vbar} ; /* mean bearing (deg) */ \ Sbar = sqrt{{Ubar*Ubar}+{Vbar*Vbar}} ; /* mean wind speed */ \ /* compute deviations from mean */ \ dU = U - Ubar ; /* */ \ dV = V - Vbar ; /* */ \ dW = W - Wbar ; /* */ \ dT = T - Tbar ; /* */ \ dA = A - Abar ; /* */ \ /* compute variances */ \ U2 = mean{dU * dU} ; /* */ \ V2 = mean{dV * dV} ; /* */ \ W2 = mean{dW * dW} ; /* */ \ T2 = mean{dT * dT} ; /* */ \ A2 = mean{dA * dA} ; /* */ \ /* compute standard deviations */ \ U1 = sqrt{ U2 } ; /* */ \ V1 = sqrt{ V2 } ; /* */ \ W1 = sqrt{ W2 } ; /* */ \ T1 = sqrt{ T2 } ; /* */ \ A1 = sqrt{ A2 } ; /* */ \ TKE = {0.5}*{sqrt{ U2+{V2+W2} }} ; /* turbulent kinetic energy */ \ /* compute covariances */ \ U1V1 = mean{ dU * dV } ; /* */ \ U1W1 = mean{ dU * dW } ; /* */ \ V1W1 = mean{ dV * dW } ; /* */ \ W1U2 = mean{ dW * {dU*dU} } ; /* */ \ W1V2 = mean{ dW * {dV*dV} } ; /* */ \ U3 = mean{ dU * {dU*dU} } ; /* */ \ V3 = mean{ dV * {dV*dV} } ; /* */ \ W3 = mean{ dW * {dW*dW} } ; /* */ \ U4 = mean{ {dU*dU} * {dU*dU} } ; /* */ \ V4 = mean{ {dV*dV} * {dV*dV} } ; /* */ \ W4 = mean{ {dW*dW} * {dW*dW} } ; /* */ \ /* set fx1 = 1 if X1 > 0.0 */ \ fU1 = requ U1 rk0 ; fU1 = {fU1-ik1}*{fU1-ik1} ; /* */ \ fV1 = requ V1 rk0 ; fV1 = {fV1-ik1}*{fV1-ik1} ; /* */ \ fW1 = requ W1 rk0 ; fW1 = {fW1-ik1}*{fW1-ik1} ; /* */ \ /* zero SKX and KUX */ \ SKU = 0.0 ; SKV = 0.0 ; SKW = 0.0 ; /* */ \ KRU = 0.0 ; KRV = 0.0 ; KRW = 0.0 ; /* */ \ if fU1 < /* compute if U1 > 0 */ \ SKU = U3 / {U1*{U1*U1}} ; /* skewness */ \ KRU = U4 / {{U1*U1}*{U1*U1}} ; /* kurtosis */ \ > ; /* */ \ if fV1 < /* compute if V1 > 0 */ \ SKV = V3 / {V1*{V1*V1}} ; /* skewness */ \ KRV = V4 / {{V1*V1}*{V1*V1}} ; /* kurtosis */ \ > ; /* */ \ if fW1 < /* compute if W1 > 0 */ \ SKW = W3 / {W1*{W1*W1}} ; /* skewness */ \ KRW = W4 / {{W1*W1}*{W1*W1}} ; /* kurtosis */ \ > ; /* */ \ /* compute temperature fluxes */ \ U1T1 = mean { dU * dT } ; /* */ \ V1T1 = mean { dV * dT } ; /* */ \ W1T1 = mean { dW * dT } ; /* */ \ /* compute buoyant prod. terms */ \ GU1T = 0.0; GV1T = 0.0; GW1T = 0.0; /* zero GX1T terms */ \ fT1 = requ Tbar rk0 ; fT1 = {fT1-ik1}*{fT1-ik1} ; /* */ \ if fT1 < /* compute if Tbar > 0 */ \ GU1T = {GRAVC * U1T1} / Tbar ; /* */ \ GV1T = {GRAVC * V1T1} / Tbar ; /* */ \ GW1T = {GRAVC * W1T1} / Tbar ; /* */ \ > ; /* */ \ /* auxilary quantities */ \ USTAR = sqrt{{U1W1*U1W1}+{V1W1*V1W1}} ; /* friction velocity */ \ /* set fT1 = 1 if W1T1 > 0.0 */ \ WSTAR = 0.0 ; MOL = 0.0 ; /* */ \ fT1 = requ W1T1 rk0 ; fT1 = {fT1-ik1}*{fT1-ik1} ; /* */ \ if fT1 < /* compute if W1T1 > 0 */ \ WSTAR = pow {sheight * GW1T} {0.3333} ; /* convection velocity scale */ \ MOL = {-1.0}*{{{USTAR*USTAR}*{USTAR*Tbar}}/{{VONKC*GRAVC}*W1T1}}; /* */ \ > ; /* */ \ SHF = DENAIR * {SHCAIR * W1T1} ; /* sensible heat flux */ \ /* save the stats for later */ \ replace Ubar tstat {vector {0} isonic} ; /* 0 - Ubar */ \ replace Vbar tstat {vector {1} isonic} ; /* 1 - Vbar */ \ replace Wbar tstat {vector {2} isonic} ; /* 2 - Wbar */ \ replace Tbar tstat {vector {3} isonic} ; /* 3 - Tbar */ \ replace Abar tstat {vector {4} isonic} ; /* 4 - Abar */ \ replace Sbar tstat {vector {5} isonic} ; /* 5 - Sbar */ \ replace U1 tstat {vector {6} isonic} ; /* 6 - Ustd */ \ replace V1 tstat {vector {7} isonic} ; /* 7 - Vstd */ \ replace W1 tstat {vector {8} isonic} ; /* 8 - Wstd */ \ replace T1 tstat {vector {9} isonic} ; /* 9 - Tstd */ \ replace A1 tstat {vector {10} isonic} ; /* 10 - Astd */ \ replace U2 tstat {vector {11} isonic} ; /* 11 - Uvar */ \ replace V2 tstat {vector {12} isonic} ; /* 12 - Vvar */ \ replace W2 tstat {vector {13} isonic} ; /* 13 - Wvar */ \ replace T2 tstat {vector {14} isonic} ; /* 14 - Tvar */ \ replace A2 tstat {vector {15} isonic} ; /* 15 - Avar */ \ replace TKE tstat {vector {16} isonic} ; /* 16 - TKE */ \ /* */ \ replace U1V1 tstat {vector {17} isonic} ; /* 17 - < U V > */ \ replace U1W1 tstat {vector {18} isonic} ; /* 18 - < U W > */ \ replace V1W1 tstat {vector {19} isonic} ; /* 19 - < V W > */ \ replace W1U2 tstat {vector {20} isonic} ; /* 20 - < W U^2 > */ \ replace W1V2 tstat {vector {21} isonic} ; /* 21 - < W V^2 > */ \ replace U3 tstat {vector {22} isonic} ; /* 22 - < U^3 > */ \ replace V3 tstat {vector {23} isonic} ; /* 23 - < V^3 > */ \ replace W3 tstat {vector {24} isonic} ; /* 24 - < W^3 > */ \ replace U4 tstat {vector {25} isonic} ; /* 25 - < U^4 > */ \ replace V4 tstat {vector {26} isonic} ; /* 26 - < V^4 > */ \ replace W4 tstat {vector {27} isonic} ; /* 27 - < W^4 > */ \ /* */ \ replace SKU tstat {vector {28} isonic} ; /* 28 - U skewness */ \ replace SKV tstat {vector {29} isonic} ; /* 29 - V skewness */ \ replace SKW tstat {vector {30} isonic} ; /* 30 - W skewness */ \ /* */ \ replace KRU tstat {vector {31} isonic} ; /* 31 - U kurtosis */ \ replace KRV tstat {vector {32} isonic} ; /* 32 - V kurtosis */ \ replace KRW tstat {vector {33} isonic} ; /* 33 - W kurtosis */ \ /* */ \ replace U1T1 tstat {vector {34} isonic} ; /* 34 - U temp flux */ \ replace V1T1 tstat {vector {35} isonic} ; /* 35 - V temp flux */ \ replace W1T1 tstat {vector {36} isonic} ; /* 36 - W temp flux */ \ /* */ \ replace GU1T tstat {vector {37} isonic} ; /* 37 - U buoyant term */ \ replace GV1T tstat {vector {38} isonic} ; /* 38 - V buoyant term */ \ replace GW1T tstat {vector {39} isonic} ; /* 39 - W buoyant term */ \ /* */ \ replace USTAR tstat {vector {40} isonic} ; /* 40 - friction veloc. */ \ replace WSTAR tstat {vector {41} isonic} ; /* 41 - conv vel scale */ \ replace MOL tstat {vector {42} isonic} ; /* 42 - Monin-Ob length */ \ replace SHF tstat {vector {43} isonic} ; /* 43 - Sen. heat flux */ \ /* */ \ nsonic = nsonic - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next sonic */ \ /* dump the results to file */ \ _dformat = "%11.5f " ; \ _dformat = "%+11.4e " ; \ N = NSONIC ; \ m = "mean(U) " ; t = extract tstat {0} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "mean{V) " ; t = extract tstat {1} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "mean(W) " ; t = extract tstat {2} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "mean(T) " ; t = extract tstat {3} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "mean(A) " ; t = extract tstat {4} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "Speed " ; t = extract tstat {5} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {6} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {7} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {8} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {9} ; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {10}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {11}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {12}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {13}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {14}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {15}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "TKE " ; t = extract tstat {16}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {17}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {18}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {19}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {20}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {21}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {22}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {23}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {24}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {25}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {26}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {27}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "SkU " ; t = extract tstat {28}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "SkV " ; t = extract tstat {29}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "SkW " ; t = extract tstat {30}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "KrU " ; t = extract tstat {31}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "KrV " ; t = extract tstat {32}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "KrW " ; t = extract tstat {33}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {34}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {35}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = " " ; t = extract tstat {36}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "g/T " ; t = extract tstat {37}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "g/T " ; t = extract tstat {38}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "g/T " ; t = extract tstat {39}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "u* " ; t = extract tstat {40}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "w* " ; t = extract tstat {41}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "L " ; t = extract tstat {42}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ m = "H " ; t = extract tstat {43}; n = N ; while n < i = N-n ; x = extract t i ; m = m + {ascii x} ; n=n-{1}; > ; print m ; \ /* */ \ report ; /* close the report file */ \ panel pid ; /* redraw the panel */ /************************************************************************* compastat - Compute the basic stats for all trials/periods *************************************************************************/ macro compastat ; \ NTRIALP = row PDATES ; /* number of trials/periods */ \ ntrialp = NTRIALP ; while ntrialp < /* loop over trials/periods */ \ itrialp = NTRIALP - ntrialp ; /* trial/period index (0:?) */ \ PCODE = itrialp ; /* set the index in the menu */ \ compbstat ; /* process the specified trial */ \ ntrialp = ntrialp - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next trial */ \ panel pid ; /* redraw the panel */ /******************************************************************************* plotcfile - Plot a corrected sonic data file Inputs: SCODE index (0:?) of the sonic file to load *************************************************************************/ macro plotcfile ; \ SFILE = extract SFILES SCODE ; /* base name of the data file */ \ SDATE = extract SDATES SCODE ; /* start date/time (y d hhmmss) */ \ SJDAY = extract SDATE {1} ; /* julian day number */ \ rfile = DSTBD + {SFILE + {".cor"}}; /* full name for corrected data */ \ pfile = DSTPD + {SFILE + {"_"}}; /* base name for plot files */ \ echo {" "} ; /* output a blank line */ \ echo {{"reading file "}+rfile} ; /* send the user a message */ \ sdata = bread rfile ; /* read the binary file */ \ sdate = mktime SDATE ; /* convert date to cal format */ \ srate = DSRATE ; /* sampling rate (Hz) */ \ nsamp = col sdata ; /* number of samples/channel */ \ tv1 = {1.0} * sdate ; /* start time (calander format) */ \ tv2 = tv1 + {{nsamp-{1.0}}/srate}; /* end time (calander format) */ \ tv = vector tv1 tv2 ; /* combine start/end time */ \ _xadate = 2 ; /* label x-axis using date/time */ \ _xlabel = "Time" ; /* define x-axis label */ \ _ymin = 0.0 ; _ymax = 1.0 ; /* set limits for null traces */ \ /* page 1 - NW(A) SW(B) Towers */ \ setwin w1 ; select {1 1} ; /* select page 1 */ \ _style = style8 ; /* draw 8-channel plot */ \ yl = "A1 - 2m AGL\nSpeed (m/s)\n" ; /* trace 1 label */ \ yl = yl + {"A1 - 2m AGL\nBearing (deg)\n"} ; /* trace 2 label */ \ yl = yl + {"A2 - 6m AGL\nSpeed (m/s)\n"} ; /* trace 3 label */ \ yl = yl + {"A2 - 6m AGL\nBearing (deg)\n"} ; /* trace 4 label */ \ yl = yl + {"B1 - 2m AGL\nSpeed (m/s)\n"} ; /* trace 5 label */ \ yl = yl + {"B1 - 2m AGL\nBearing (deg)\n"} ; /* trace 6 label */ \ yl = yl + {"B2 - 6m AGL\nSpeed (m/s)\n"} ; /* trace 7 label */ \ yl = yl + {"B2 - 6m AGL\nBearing (deg)"} ; /* trace 8 label */ \ _ylabel = yl ; /* y-axis label */ \ _title = "DSTL Sonics - West Towers - Day " ; /* define the title */ \ _title = _title + {{ascii SJDAY} + {" - "}}; /* add the julian day */ \ _title = _title + {asctime sdate} ; /* add the start time */ \ UA1 = block {extract sdata { 0}} NBLOCK ; /* A1 - U component (m/s) */ \ VA1 = block {extract sdata { 1}} NBLOCK ; /* A1 - V component (m/s) */ \ UA2 = block {extract sdata { 4}} NBLOCK ; /* A2 - U component (m/s) */ \ VA2 = block {extract sdata { 5}} NBLOCK ; /* A2 - V component (m/s) */ \ UB1 = block {extract sdata { 8}} NBLOCK ; /* B1 - U component (m/s) */ \ VB1 = block {extract sdata { 9}} NBLOCK ; /* B1 - V component (m/s) */ \ UB2 = block {extract sdata {12}} NBLOCK ; /* B2 - U component (m/s) */ \ VB2 = block {extract sdata {13}} NBLOCK ; /* B2 - V component (m/s) */ \ SA1 = sqrt{{UA1*UA1}+{VA1*VA1}} ; /* A1 - wind speed (m/s) */ \ SA2 = sqrt{{UA2*UA2}+{VA2*VA2}} ; /* A2 - wind speed (m/s) */ \ SB1 = sqrt{{UB1*UB1}+{VB1*VB1}} ; /* B1 - wind speed (m/s) */ \ SB2 = sqrt{{UB2*UB2}+{VB2*VB2}} ; /* B2 - wind speed (m/s) */ \ BA1 = RADTOD * {atan UA1 VA1} ; /* A1 - wind bearing (deg) */ \ BA2 = RADTOD * {atan UA2 VA2} ; /* A2 - wind bearing (deg) */ \ BB1 = RADTOD * {atan UB1 VB1} ; /* B1 - wind bearing (deg) */ \ BB2 = RADTOD * {atan UB2 VB2} ; /* B2 - wind bearing (deg) */ \ y = matrix SA1 BA1 SA2 BA2 SB1 BB1 SB2 BB2; /* combine wind vectors */ \ gid1 = draw tv y ; /* draw the series */ \ setwin w1 ; select {1 2} ; /* select page 2 */ \ y = matrix UA1 VA1 UA2 VA2 UB1 VB1 UB2 VB2; /* combine wind vectors */ \ gid2 = draw tv y ; /* draw the series */ \ /* page 3 - NE(C) SE(D) Towers */ \ setwin w1 ; select {1 3} ; /* select page 3 */ \ _style = style8 ; /* draw 8-channel plot */ \ yl = "C1 - 2m AGL\nSpeed (m/s)\n" ; /* trace 1 label */ \ yl = yl + {"C1 - 2m AGL\nBearing (deg)\n"} ; /* trace 2 label */ \ yl = yl + {"C2 - 6m AGL\nSpeed (m/s)\n"} ; /* trace 3 label */ \ yl = yl + {"C2 - 6m AGL\nBearing (deg)\n"} ; /* trace 4 label */ \ yl = yl + {"D1 - 2m AGL\nSpeed (m/s)\n"} ; /* trace 5 label */ \ yl = yl + {"D1 - 2m AGL\nBearing (deg)\n"} ; /* trace 6 label */ \ yl = yl + {"D2 - 6m AGL\nSpeed (m/s)\n"} ; /* trace 7 label */ \ yl = yl + {"D2 - 6m AGL\nBearing (deg)"} ; /* trace 8 label */ \ _ylabel = yl ; /* y-axis label */ \ _title = "DSTL Sonics - East Towers - Day " ; /* define the title */ \ _title = _title + {{ascii SJDAY} + {" - "}}; /* add the julian day */ \ _title = _title + {asctime sdate} ; /* add the start time */ \ UC1 = block {extract sdata {16}} NBLOCK ; /* C1 - U component (m/s) */ \ VC1 = block {extract sdata {17}} NBLOCK ; /* C1 - V component (m/s) */ \ UC2 = block {extract sdata {20}} NBLOCK ; /* C2 - U component (m/s) */ \ VC2 = block {extract sdata {21}} NBLOCK ; /* C2 - V component (m/s) */ \ UD1 = block {extract sdata {24}} NBLOCK ; /* D1 - U component (m/s) */ \ VD1 = block {extract sdata {25}} NBLOCK ; /* D1 - V component (m/s) */ \ UD2 = block {extract sdata {28}} NBLOCK ; /* D2 - U component (m/s) */ \ VD2 = block {extract sdata {29}} NBLOCK ; /* D2 - V component (m/s) */ \ SC1 = sqrt{{UC1*UC1}+{VC1*VC1}} ; /* C1 - wind speed (m/s) */ \ SC2 = sqrt{{UC2*UC2}+{VC2*VC2}} ; /* C2 - wind speed (m/s) */ \ SD1 = sqrt{{UD1*UD1}+{VD1*VD1}} ; /* D1 - wind speed (m/s) */ \ SD2 = sqrt{{UD2*UD2}+{VD2*VD2}} ; /* D2 - wind speed (m/s) */ \ BC1 = RADTOD * {atan UC1 VC1} ; /* C1 - wind bearing (deg) */ \ BC2 = RADTOD * {atan UC2 VC2} ; /* C2 - wind bearing (deg) */ \ BD1 = RADTOD * {atan UD1 VD1} ; /* D1 - wind bearing (deg) */ \ BD2 = RADTOD * {atan UD2 VD2} ; /* D2 - wind bearing (deg) */ \ y = matrix SC1 BC1 SC2 BC2 SD1 BD1 SD2 BD2; /* combine wind vectors */ \ gid2 = draw tv y ; /* draw the series */ \ setwin w1 ; select {1 4} ; /* select page 4 */ \ y = matrix UC1 VC1 UC2 VC2 UD1 VD1 UD2 VD2; /* combine wind vectors */ \ gid4 = draw tv y ; /* draw the series */ \ /* create Postscript plots */ \ _pfile = pfile + {"W.ps"} ; plot {1 1} ; /* West towers */ \ _pfile = pfile + {"E.ps"} ; plot {1 3} ; /* East towers */ \ panel pid ; /* redraw the panel */ /******************************************************************************* makeafile - Create ascii files from a corrected sonic data file Inputs: SCODE index (0:?) of the sonic file to load *************************************************************************/ macro makeofile ; \ SFILE = extract SFILES SCODE ; /* base name of the data file */ \ SDATE = extract SDATES SCODE ; /* start date/time (y d hhmmss) */ \ SJDAY = extract SDATE {1} ; /* julian day number */ \ rfile = DSTBD + {SFILE + {".cor"}}; /* full name for corrected data */ \ afile = DSTAD + {SFILE + {"_"}}; /* base name for ascii files */ \ echo {" "} ; /* output a blank line */ \ echo {{"reading file "}+rfile} ; /* send the user a message */ \ sdata = bread rfile ; /* read the binary file */ \ sdate = mktime SDATE ; /* convert date to cal format */ \ srate = DSRATE ; /* sampling rate (Hz) */ \ nsamp = col sdata ; /* number of samples/channel */ \ rv1 = 0.0 ; /* relative start time */ \ rv2 = rv1 + {{nsamp-{1.0}}/srate}; /* end time */ \ rv = line rv1 rv2 nsamp ; /* time vector */ \ _ncol = 5 ; /* output 5 column (time, U,..) */ \ _dformat = "%.3f" ; /* output format */ \ nsonic = NSONIC ; while nsonic < /* loop over sonics */ \ isonic = NSONIC - nsonic ; /* sonic index (0:NSONIC-1) */ \ soname = extract SNAMES isonic ; /* name of the sonic (N1,...) */ \ ofile = afile + {soname+{".dat"}}; /* full name of the output file */ \ echo {{"creating file "}+ofile} ; /* send the user a message */ \ r1 = {4}*isonic ; /* first row to extract */ \ r2 = r1 + {3} ; /* last row to extract */ \ c1 = 0 ; /* first col to extract */ \ c2 = nsamp - {1} ; /* last col to extract */ \ rc1 = vector r1 c1 ; /* combine first row/col */ \ rc2 = vector r2 c2 ; /* combine last row/col */ \ tdata = extract sdata rc1 rc2 ; /* extract current sonic data */ \ U = extract tdata {0}; /* x component (m/s) */ \ V = extract tdata {1}; /* y component (m/s) */ \ W = extract tdata {2}; /* z component (m/s) */ \ T = extract tdata {3}; /* temperature (K) */ \ odata = matrix rv U V W T ; /* form output matrix */ \ odata = transpose odata ; /* convert to 5-column matrix */ \ awrite odata ofile ; /* write the matrix to file */ \ nsonic = nsonic - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next sonic */ \ panel pid ; /* redraw the panel */ macro makeafile ; \ NSFILEP = row SFILES ; /* number of sonic files */ \ nsonicp = NSFILEP ; while nsonicp < /* loop over sonic files */ \ isonicp = NSFILEP - nsonicp ; /* sonics file index (0:?) */ \ SCODE = isonicp ; /* set the index in the menu */ \ makeofile ; /* process the specified file */ \ nsonicp = nsonicp - {1} ; /* decrement the loop counter */ \ > ; /* continue with the next file */ \ panel pid ; /* redraw the panel */ /************************************************************************* Part 3 - switch to graphics mode, ... *************************************************************************/ _tsize = TWINDSZ /* origin/size of the terminal */ TEMPLATE = "dstlsonic.pan" /* name of the menu file */ switch /* switch to graphics mode */ setpal rgb /* set the color palette */ w0= openwin {"SONIC PROCESSOR"} PWINDSZ /* open the menu */ w1= openwin {"SONIC DATA"} GWINDSZ {1} /* open the display */ setwin w0 ; select {1 1} /* select page 1 of w0 */ pid = panel TEMPLATE ; show {1 1} /* load and display the panel */ setwin w1 ; select {1 1} /* select page 1 of w1 */ free ASPECT GWINDSZ LMARGIN PGHEIGHT /* release temporary variables */ free PGWIDTH PHEIGHT PWIDTH PORIENT /* release temporary variables */ free PRRES PRTYPE PWINDSZ SHOWPAGE /* release temporary variables */ free TEMPLATE TMARGIN TWINDSZ /* release temporary variables */ /* wait for user input */