OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
read_gps.f
Go to the documentation of this file.
1  subroutine read_gps(input,nframes,scal_p,xmaglm,gpsvec,nsig,
2  * igyr,igday,gpsec,secst,secend,ngps)
3 c
4 c Purpose: This subroutine extracts the valid GPS data from the input
5 c telemetry structure. Only data points with four or more
6 c GPS signals are included. Also, the time tags for all
7 c input minor frames are checked to determine the time
8 c range of the data.
9 c
10 c Calling Arguments:
11 c
12 c Name Type I/O Description
13 c -------- ---- --- -----------
14 c input struct I input data structure
15 c nframes I*4 I Number of frames in input structure
16 c scal_p R*8 I Scale factor for position vectors
17 c xmaglm(2) R*8 I Limits on position magnitude
18 c gpsvec(6,*) R*8 O Output GPS orbit vectors
19 c nsig(*) I*4 O Number of GPS signals for each vector
20 c igyr I*4 O Time tag year
21 c igday I*4 O Time tag day-of-year
22 c gpsec(*) R*8 O Time tag seconds of igday
23 c secend R*8 O End of interval time in seconds
24 c ngps I*4 O Number of GPS vectors
25 c
26 c By: Frederick S. Patt, GSC, December 21, 1993
27 c
28 c Notes:
29 c
30 c Modification History:
31 c
32 c Modified to unpack GPS time tag from telemetry (previously used
33 c minor frame time tag). F. S. Patt, GSC, September 14, 1996
34 c
35 c Modified to perform sanity check on time tags (e.g., all tags
36 c fall within +/- 1 day of first valid time tag.
37 c F. S. Patt, GSC, October 4, 1996.
38 c
39 c Fixed a bug in validating the first GPS time code.
40 c F. S. Patt, SAIC GSC, September 7, 1997.
41 c
42 c Added additional range checking to the seconds-of-day.
43 c F. S. Patt, SAIC GSC, September 10, 1997
44 c
45 c Changed integer*1 to byte for Sun OS compatibility, B. A. Franz,
46 c GSC, November 14, 1997.
47 c
48 c Modified to reject GPS samples with time tags later than minor frame.
49 c F. S. Patt, SAIC GSC, July 22, 1998.
50 c
51 c Modified to reject GPS samples prior to current time range.
52 c F. S. Patt, SAIC GSC, March 3, 2000
53 c
54 c Modified to remove checks on minor frame times (SWl01 takes care of this).
55 c F. S. Patt, SAIC GSC, May 20, 2000
56 c
57 c Fixed bug which caused error message to be output for duplicate GPS sample.
58 c F. S. Patt, SAIC GSC, June 8, 2000
59 c
60 c Fixed bug introduced in indexing of input array during generation of
61 c single-source code for SGI and Linux.
62 c F. S. Patt, SAIC GSC, January 19, 2001
63 c
64 c Modified duplicate time check to use a tolerance instead of equivalence.
65 c F. S. Patt, SAIC, February 17, 2004
66 
67  implicit none
68 #include "nav_cnst.fin"
69 #include "input_s.fin"
70 
71  type(input_struct) :: input(maxlin)
72 
73  real*8 gpsvec(6,maxlin), gpsec(maxlin), scal_p, secst, secend
74  real*8 sec, xmaglm(2), pmag
75  integer*4 nsig(maxlin), nframes, igyr, igday, ngps
76  integer*4 iframe, mframe, if1, j, iy2, im2, id2, ms, mp
77  integer*4 jd, jd1, jd2, jdtol
78  integer*2 ib2
79  byte b2(2)
80  logical gotone,gotime
81  data jdtol/1/
82  equivalence(ib2,b2)
83 
84 c Find two good minor frames with the same day in the time tag
85  if1 = 1
86  ngps = 0
87  gotime = .true.
88  igyr = input(if1)%iyear
89  igday = input(if1)%iday
90  jd1 = jd(igyr,1,igday)
91  secst = input(if1)%msec/1000.d0
92  secend = secst
93 
94 c Find first good vector
95  gotone = .false.
96  dowhile(.not.gotone)
97  mframe = input(if1)%sc_id(1)/128
98 
99 c Check minor frame time tag against time range
100  iy2 = input(if1)%iyear
101  id2 = input(if1)%iday
102  jd2 = jd(iy2,1,id2)
103  sec = input(if1)%msec/1000.d0
104  sec = sec +(jd2-jd1)*864.d2
105  if (sec.lt.secst) secst = sec
106  if (sec.gt.secend) secend = sec
107  if((input(if1)%flag.eq.0).and.(mframe.eq.1)) then
108 
109 c Unpack number of GPS signals from telemetry data and check for minimum
110  ms = input(if1)%sc_dis(12)
111  if (ms.ge.4) then
112 
113 c Get GPS time tag
114 c Year is two bytes
115 #ifdef LINUX
116  b2(1) = input(if1)%sc_dis(15)
117  b2(2) = input(if1)%sc_dis(14)
118 #else
119  b2(1) = input(if1)%sc_dis(14)
120  b2(2) = input(if1)%sc_dis(15)
121 #endif
122  iy2 = ib2
123 c Month and day
124  im2 = input(if1)%sc_dis(16)
125  id2 = input(if1)%sc_dis(17)
126  jd2 = jd(iy2,im2,id2)
127 c Convert hours, minutes, seconds, fraction to seconds of day
128  sec = input(if1)%sc_dis(18)*36.d2
129  * + input(if1)%sc_dis(19)*6.d1 + input(if1)%sc_dis(20)
130  * + input(if1)%sc_ana(21)
131 
132  sec = sec + (jd2-jd1)*864.d2
133 
134 
135  if ( (sec.ge.(secst-10.)) .and. (sec.lt.secend) ) then
136 
137  ngps = ngps + 1
138  gotone = .true.
139  if (sec.lt.secst) secst = sec
140 
141 c Store data in output array
142  gpsec(ngps) = sec
143  nsig(ngps) = ms
144  do j=1,3
145  gpsvec(j,ngps) = input(if1)%sc_ana(j)*scal_p
146  gpsvec(j+3,ngps) = input(if1)%sc_ana(j+3)
147  end do
148 c Perform limit check on position magnitude
149  pmag = sqrt(gpsvec(1,ngps)**2 + gpsvec(2,ngps)**2
150  * + gpsvec(3,ngps)**2)
151  if ((pmag.lt.xmaglm(1)).or.(pmag.gt.xmaglm(2))) then
152  gotone = .false.
153  ngps = ngps - 1
154  end if
155  end if
156  end if
157  end if
158  if1 = if1 + 1
159  if (if1.gt.nframes) go to 999
160  end do
161  mp = ms
162 
163 c Now unpack remaining GPS vectors; check for >4 signals at previous point.
164  do iframe = if1,nframes
165  iy2 = input(iframe)%iyear
166  id2 = input(iframe)%iday
167  jd2 = jd(iy2,1,id2)
168  sec = input(iframe)%msec/1000.d0
169  sec = sec + (jd2-jd1)*864.d2
170 
171 c Check time tag against range
172  if (sec.lt.secst) secst = sec
173  if (sec.gt.secend) secend = sec
174 
175 c Check for minor frame 1 and valid data
176  mframe = input(iframe)%sc_id(1)/128
177  if ((input(iframe)%flag.eq.0).and.(mframe.eq.1)) then
178 
179 c Unpack number of GPS signals from telemetry data
180  ms = input(iframe)%sc_dis(12)
181 
182 c Get GPS time tag
183 c Year is two bytes
184 #ifdef LINUX
185  b2(1) = input(iframe)%sc_dis(15)
186  b2(2) = input(iframe)%sc_dis(14)
187 #else
188  b2(1) = input(iframe)%sc_dis(14)
189  b2(2) = input(iframe)%sc_dis(15)
190 #endif
191  iy2 = ib2
192 c Month and day
193  im2 = input(iframe)%sc_dis(16)
194  id2 = input(iframe)%sc_dis(17)
195  jd2 = jd(iy2,im2,id2)
196  sec = input(iframe)%sc_dis(18)*36.d2
197  * + input(iframe)%sc_dis(19)*6.d1
198  * + input(iframe)%sc_dis(20)
199  * + input(iframe)%sc_ana(21)
200 
201 c Convert hours, minutes, seconds, fraction to seconds of day
202  sec = sec + (jd2-jd1)*864.d2
203 
204 c Check for time tag out of range
205  if ((sec.ge.secst).and.(sec.le.secend)) then
206 
207 c Check for < minimum GPS signals or duplicated time tag
208  if ((ms.ge.4).and.(mp.ge.4).and.
209  * (dabs(sec-gpsec(ngps)).gt.1.d-9)) then
210  ngps = ngps + 1
211 
212 c Store data in output array
213  gpsec(ngps) = sec
214  nsig(ngps) = ms
215  do j=1,3
216  gpsvec(j,ngps) = input(iframe)%sc_ana(j)*scal_p
217  gpsvec(j+3,ngps) = input(iframe)%sc_ana(j+3)
218  end do
219 c Perform limit check on position magnitude
220  pmag = sqrt(gpsvec(1,ngps)**2 + gpsvec(2,ngps)**2
221  * + gpsvec(3,ngps)**2)
222  if ((pmag.lt.xmaglm(1)).or.(pmag.gt.xmaglm(2))) then
223  ngps = ngps - 1
224  end if
225 
226 c Save number of GPS signals
227  mp = ms
228  end if
229  else
230  print *,'READ_GPS: Invalid date',iy2,im2,id2,
231  * ' record',iframe
232  end if
233  end if
234  end do
235 
236  999 print *,'READ_GPS:',ngps,' valid GPS vectors'
237  print *,' starting at ',igyr,igday,secst
238  print *,' ending at ',igyr,igday,secend
239 
240  return
241  end
#define real
Definition: DbAlgOcean.cpp:26
subroutine read_gps(input, nframes, scal_p, xmaglm, gpsvec, nsig, igyr, igday, gpsec, secst, secend, ngps)
Definition: read_gps.f:3
Definition: jd.py:1