program findvectors c - 2018 03 21 c - Originally created for NADCON 5, this version has been modified c - to work for VERTCON 3 c - Updated 2015 10 07 for new vector format c - Analytical program built during the NADCON 5.0 process c - Built to scan a GMT-based vector file and either: c - 1) Spit out all vectors near the input lat/lon c - 2) Spit out *outlier* vectors near the input lat/lon c - Expected format would look like this: c c - VERTCON 3: c LON LAT Meters PID c (rvalm) c 278.2166665250 24.5500001417 -0.415 AA0001 c 278.1999999639 24.5500001111 -0.408 AA0002 c 278.1946555556 24.5606000000 -0.408 AA0003 c 278.1938055556 24.5600000000 -0.409 AA0004 implicit real*8(a-h,o-z) parameter(max = 500000) character*80 card character*200 fname character*1 type real*8 glo(max),gla(max),az(max),sval(max) real*8 rvalm(max) character*6 xpid,pid(max) logical*1 ldrop(max),noaz write(6,1) 1 format('Program findvectors') write(6,2) 2 format('Enter file name to be searched: ',$) read(5,'(a)')fname write(6,3) 3 format('Enter central lat/lon for search : ',$) read(5,*)xlat0,xlon0 write(6,4) 4 format('Enter search distance in arcdegrees : ',$) read(5,*)dx 6 write(6,5) 5 format('Searching for (a)ll or (o)utliers? : ',$) read(5,'(a1)')type if(type.ne.'a' .and. type.ne.'o')goto 6 if(type.eq.'o')then write(6,10) 10 format(' Enter outlier sigma multiplier: ',$) read(5,*)SigmaMultiplier endif write(6,7)trim(fname) 7 format('Opening file: ',a) open(1,file=fname,status='old',form='formatted') i = 0 j = 0 c 101 format(f16.10,1x,f15.10,1x,f6.2,1x, c *f12.2,1x,f9.5,1x,f9.3,1x,a6) 101 format(f16.10,1x,f15.10,1x,6x ,1x, *12x ,1x,9x ,1x,f9.3,1x,a6) c 103 format('Bad Vector Length ', c *f16.10,1x,f15.10,1x,f6.2,1x, c *f12.2,1x,f9.5,1x,f9.3,1x,a6) 103 format('Bad Vector Length ', *f16.10,1x,f15.10,1x,6x ,1x, *12x ,1x,9x ,1x,f9.3,1x,a6) c 100 read(1,101,end=102)xlo,xla,xaz,xsval,xrvals,xrvalm,xpid 100 read(1,101,end=102)xlo,xla ,xrvalm,xpid i = i + 1 dla = dabs(xla - xlat0) if(dla.gt.dx)goto 100 dlo = dabs(xlo - xlon0) if(dlo.gt.dx)goto 100 j = j + 1 glo(j) = xlo gla(j) = xla c az(j) = xaz c sval(j)= xsval c rvals(j)= xrvals rvalm(j)= xrvalm pid(j) = xpid ldrop(j) = .false. goto 100 102 npts = i ngpts = j write(6,8)npts,ngpts 8 format('Done reading file. ',/, *6x,'Total pts found : ',i10,/, *6x,'Total pts found in search box: ',i10) 301 format('All points in search box: ') if(type.eq.'a')then write(6,301) do 300 i=1,ngpts c write(6,101)glo(i),gla(i),az(i),sval(i), c * rvals(i),rvalm(i),pid(i) write(6,101)glo(i),gla(i), * rvalm(i),pid(i) 300 continue stop endif c - At this point, type must be "o" iter = 0 800 iter = iter + 1 write(6,801)iter 801 format('Iteration # ',i8) c aveaz = 0.d0 c rmsaz = 0.d0 c - For VERTCON 3, only working on OHT (meters) c - Do all searches on METERS, as lat/lon/eht all have them, but c - eht does not have SECONDS. averv= 0.d0 rmsrv= 0.d0 ikt = 0 do 200 i=1,ngpts if(ldrop(i))goto 200 c aveaz = aveaz + az(i) c rmsaz = rmsaz + az(i)**2 averv = averv + rvalm(i) rmsrv = rmsrv + rvalm(i)**2 ikt = ikt + 1 200 continue c aveaz = aveaz / ikt c rmsaz = dsqrt(rmsaz / ikt) averv = averv / ikt rmsrv = dsqrt(rmsrv / ikt) fact = dble(ikt) / dble(ikt-1) c stdaz = dsqrt(fact*(rmsaz**2 - aveaz**2)) stdrv = dsqrt(fact*(rmsrv**2 - averv**2)) c write(6,202)aveaz,stdaz,averv,stdrv write(6,202) averv,stdrv c 202 format('Searching for outliers deviating from:',/, c *6x,'Ave Az = ',f0.5,/, c *6x,'Std Az = ',f0.5,/, c *6x,'Ave Val(m) = ',f0.5,/, c *6x,'Std Val(m) = ',f0.5) 202 format('Searching for outliers deviating from:',/, *6x,'Ave Val(m) = ',f0.5,/, *6x,'Std Val(m) = ',f0.5) rvlo = averv - SigmaMultiplier*stdrv rvhi = averv + SigmaMultiplier*stdrv c azlo = aveaz - SigmaMultiplier*stdaz c azhi = aveaz + SigmaMultiplier*stdaz c write(6,803)rvlo,rvhi,azlo,azhi write(6,803)rvlo,rvhi c 803 format('Real Value, Low (m): ',f0.5,/, c * 'Real Value, High(m): ',f0.5,/, c * 'Azimuth , Low : ',f0.5,/, c * 'Azimuth , High : ',f0.5) 803 format('Real Value, Low (m): ',f0.5,/, * 'Real Value, High(m): ',f0.5) noaz = .false. c if(SigmaMultiplier * stdaz .gt. 180)then c write(6,11) c 11 format('Azimuth Sigma too large to ', c * 'support outlier detection. No search on Azimuth') c noaz = .true. c endif inew = 0 do 201 i=1,ngpts if(ldrop(i))goto 201 c - Drop on Real Value if(rvalm(i).lt.rvlo .or. rvalm(i).gt.rvhi)then c write(6,103)glo(i),gla(i),az(i),sval(i), c * rvals(i),rvalm(i),pid(i) write(6,103)glo(i),gla(i), * rvalm(i),pid(i) ldrop(i) = .true. inew = inew + 1 goto 201 endif if(noaz)goto 888 c - Drop on Azimuth (caution...wrap-around issue...) c - Case 1: Low azimuth drops below zero c if(azlo.lt.0)then c azlo1 = azlo + 360.d0 c if(az(i).lt.azlo1 .and. az(i).gt.azhi)then c write(6,104)glo(i),gla(i),az(i),sval(i), c * rvals(i),rvalm(i),pid(i) c ldrop(i) = .true. c inew = inew + 1 c endif c - Case 2: High azimuth rises above 360 c elseif(azhi.gt.360)then c azhi1 = azhi - 360.d0 c if(az(i).lt.azlo .and. az(i).gt.azhi1)then c write(6,104)glo(i),gla(i),az(i),sval(i), c * rvals(i),rvalm(i),pid(i) c ldrop(i) = .true. c inew = inew + 1 c endif c - Case 3: Both high and low azimuth remain between 0 and 360 c else c if(az(i).lt.azlo .or. az(i).gt.azhi)then c write(6,104)glo(i),gla(i),az(i),sval(i), c * rvals(i),rvalm(i),pid(i) c ldrop(i) = .true. c inew = inew + 1 c endif c endif 888 continue 201 continue if(inew.ne.0)then write(6,802)inew goto 800 endif 802 format('New outliers found this iteration: ',i8) end