NADCON5-ng  0.0.1
NADCON5 Next Generation
makework.f
Go to the documentation of this file.
1 c> \ingroup doers
2 c> **Program** to create a *work* file which will
3 c> serve as the primary information needed to
4 c> analyze and create NADCON v5.0 grids.
5 c>
6 c>
7 c> This program is based on previous programs that
8 c> were created by Dru Smith during the GEOCON v2.0
9 c> process. It has been modified specifically to be a tool
10 c> for NADCON v5.0.
11 c>
12 c> Rather than have multiple programs (1, 2, 3, 4 as
13 c> was the case for GEOCON v2.0), It was decided
14 c> make ONE working "makework.f" program (this one)
15 c> and have it feed off of an input file which
16 c> can be modified.
17 c>
18 c> The input file will reflect all that is necessary
19 c> to create a work file.
20 c>
21 c> ### Program arguments
22 c> Arguments are newline terminated and read from standard input
23 c>
24 c> They are enumerated here
25 c> \param oldtm Source Datum
26 c> \param newdtm Target Datum,region
27 c> \param region Conversion Region
28 c>
29 c> ### Program Inputs:
30 c>
31 c>- A *control file* in directory `Control/`, the name is generated from
32 c> input arguments
33 c>
34 c> Known control file names are:
35 c>
36 c> cfname = control.ussd.nad27.conus
37 c> cfname = control.nad27.nad83_1986.conus
38 c>
39 c>- A *manual edits* file, called `workedits` in
40 c> directory `Work/`
41 c>
42 c> ## By way of example...
43 c>
44 c> **If** the input file controlling the creation of the work file is:
45 c>
46 c> Control/control.ussd.nad27.conus
47 c>
48 c> then the output data file is:
49 c>
50 c> Work/work.ussd.nad27.conus
51 c>
52 c> The work file has the following format:
53 c>
54 c> Cols Format Description
55 c> 1- 6 a6 PID
56 c> 7 1x - blank -
57 c> 8- 9 a2 State
58 c> 10 a1 Reject code for missing latitude pair (blank for good)
59 c> 11 a1 Reject code for missing longitude pair (blank for good)
60 c> 12 a1 Reject code for missing ellip ht pair (blank for good)
61 c> 13 1x - blank -
62 c> 14- 27 f14.10 Latitude (Old Datum), decimal degrees (-90 to +90)
63 c> 28 1x - blank -
64 c> 29- 42 f14.10 Lonitude (Old Datum), decimal degrees (0 to 360)
65 c> 43 1x - blank -
66 c> 44- 51 f8.3 Ellipsoid Height (Old datum), meters
67 c> 52 1x - blank -
68 c> 53- 61 f9.5 Delta Lat (New Datum minus Old Datum), arcseconds
69 c> 62 1x - blank -
70 c> 63- 71 f9.5 Delta Lon (New Datum minus Old Datum), arcseconds
71 c> 72 1x - blank -
72 c> 73- 81 f9.3 Delta Ell Ht (New Datum minus Old Datum), meters
73 c> 82 1x - blank -
74 c> 83- 91 f9.3 Delta Horizontal (absolute value), arcseconds
75 c> 92 1x - blank -
76 c> 93-101 f9.5 Azimuth of Delta Horizontal (0-360), degrees
77 c> 102 1x - blank -
78 c> 103-111 f9.3 Delta Lat (New Datum minus Old Datum), meters
79 c> 112 1x - blank -
80 c> 113-121 f9.3 Delta Lon (New Datum minus Old Datum), meters
81 c> 122 1x - blank -
82 c> 123-131 f9.3 Delta Horizontal (absolute value), meters
83 c> 132 1x - blank -
84 c> 133-142 a10 Old Datum Name
85 c> 143 1x - blank -
86 c> 144-153 a10 New Datum Name
87 c> format(a6,1x,a2,a1,a1,a1,1x,f14.10,1x,f14.10,1x,f8.3,1x,
88 c> *f9.5,1x,f9.5,1x,f9.3,1x,f9.3,1x,f9.5,1x,f9.3,1x,f9.3,1x,f9.3,
89 c> 1x,a10,1x,a10)
90 c>
91 c>
92 c> This differs from Dennis's GEOCON v1.0 in that:
93 c> - 3 reject codes
94 c> - 10 decimal places in latitude (See DRU-10, p. 123)
95 c> - 10 decimal places in longitude (See DRU-10, p. 123)
96 c> - Lat, Lon and Horizontal in both arcseconds and meters each
97 c> - Azimuth of Horizontal
98 c> - Identification of which datums are transformed
99 c>
100 c> ##References
101 c> ### NADCON5:
102 c>
103 c> See:
104 c> - DRU-11, p. 124
105 c>
106 c> ### GEOCON v2.0:
107 c>
108 c> See:
109 c> - DRU-10, p. 143
110 c> - DRU-11, p. 10
111 c> - DRU-11, p. 26
112 c> - DRU-11, p. 56
113 c>
114 c> ## Changelog
115 c>
116 c> ### 2017 11 19 (**NG**)
117 c> Formated Comments to be compatible with Doxygen
118 c> Due to deprecation of various arguments for GMT tools, invocation of `xyz2grd`
119 c> and `grd2xyz` have arguments options changed in generated files
120 c> - `-bos` -> `-bo3f` with equivalent meaning of a 3 column single precision output
121 c> - `-bis` -> `-bi3f` with equivalent meaning of a 3 column single precision input
122 c>
123 c> ### 2016 09 14:
124 c> Due to some complications in mymedian5.f which happen if
125 c> data is in the *work* file that is not to be sorted and used, I've decided
126 c> to put the *out of grid* point removal code here, so that such points
127 c> will go into the *work* file but will all get a `111` set of reject codes
128 c> so they don't go forward in the processing.
129 c>
130 c> ### 2016 09 13:
131 c> Fixed a bug that sends an incoming `0` reject flag as a **zero**. All
132 c> later programs expect a BLANK for a "good" reject code. The incoming `0` is
133 c> from the `workedits` file and is fine to come in, but must go OUT as a BLANK.
134 c>
135 c> Also, put in code to correct situations where an entry is in `workedits`, but
136 c> the PID for that entry isn't actually in the incoming data. This relies
137 c> on a new vector `EditTracker`
138 c>
139 c> ### 2016 02 26:
140 c> Change (see: DRU-12, p. 18) to reflect the decision that "manual edits" should ONLY
141 c> edit data OUT and **not** add data back in.
142 c>
143 c> ### 2016 01 07:
144 c> Changed to split "Relevant Edits" into
145 c> three counts: lat, lon and eht
146  program makework
147 
148  implicit double precision(a-h,o-z)
149  parameter(maxedits = 10000)
150 
151  character*58 fname0
152  character*200 cfname
153  character*200 wfname
154  character*200 efname
155  character*80 card
156 
157 c - Information extracted from control file
158  character*6 cline
159  character*80 cheader
160  character*30 cregion
161  character*30 cdatum1
162  character*30 cdatum2
163 
164 c - Variables for working with manual edit file
165  character*10 dummy1,dummy2,dummy3
166  character*10 EditRegion(maxedits)
167  character*10 EditOldDtm(maxedits)
168  character*10 EditNewDtm(maxedits)
169  character*6 EditPID(maxedits)
170  character*1 EditRejLat(maxedits)
171  character*1 EditRejLon(maxedits)
172  character*1 EditRejEht(maxedits)
173 
174 c - Each *.in file (including directory):
175  character*100 fname
176 
177 c - "h" and "f" are artifacts. "h" for HARN meaning
178 c - "old datum". "f" for FBN meaning "new datum".
179 c - Too much trouble to go back and fix variable names.
180 c - Just know that "f" means NEW and "h" means OLD,
181 c - no matter what the datums themselves ARE.
182  character*15 nameh
183  character*15 namef
184  character*6 pid
185  character*2 state
186  character*13 clath,clatf
187  character*14 clonh,clonf
188  character*9 cehth,cehtf
189 
190  character*10 olddtm,newdtm,region
191 
192  character*1 rejlat,rejlon,rejeht
193 
194  character*200 suffix1
195 
196  logical badlat,badlon,badeht
197 
198 c - 2016 09 13 - To track the use of each so-called "relevant" edit,
199 c - in case some of them aren't actually relevant (e.g. the PID
200 c - doesn't match what's coming in)
201  logical EditTracker(maxedits)
202 
203 c ------------------------------------------------------------------
204 c - BEGIN PROGRAM
205 c ------------------------------------------------------------------
206  write(6,1001)
207  1001 format('BEGIN program makework.f')
208 
209 c ------------------------------------------------------------------
210 c - User-supplied input
211 c ------------------------------------------------------------------
212  read(5,'(a)')olddtm
213  read(5,'(a)')newdtm
214  read(5,'(a)')region
215 
216 c ------------------------------------------------------------------
217 c - 2016 09 14 : Get the official grid bounds for this region,
218 c - so that we can auto-reject (Flag with "111") any points
219 c - that come in and are not inside the official grid bounds.
220 c ------------------------------------------------------------------
221  call getgridbounds(region,xn,xs,xw,xe)
222  noutside = 0
223 
224 c ------------------------------------------------------------------
225 c - Generate the suffixes used in all our files
226 c ------------------------------------------------------------------
227  suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region)
228 
229 c ------------------------------------------------------------------
230 c - Open the control file
231 c ------------------------------------------------------------------
232  cfname='Control/control.'//trim(suffix1)
233  open(1,file=cfname,status='old',form='formatted')
234  write(6,1004)trim(cfname)
235  1004 format(6x,'makework.f: Accessing control file ',a)
236 
237 c ------------------------------------------------------------------
238 c - Open the "manual edits" file
239 c ------------------------------------------------------------------
240  efname='Work/workedits'
241  open(20,file=efname,status='old',form='formatted')
242  write(6,1006)trim(efname)
243  1006 format(6x,'makework.f: Accessing workedits file ',a)
244 
245 c ------------------------------------------------------------------
246 c - Create and open the work file
247 c ------------------------------------------------------------------
248 c - The location of the "work*" file
249  wfname = 'Work/work.'//trim(suffix1)
250  open(2,file=wfname,status='new',form='formatted')
251 
252  write(6,1002)trim(wfname)
253  1002 format(6x,'makework.f: Creating work file ',a)
254 
255 c ------------------------------------------------------------------
256 c - Some necessary constants.
257 c ------------------------------------------------------------------
258  pi = 2.d0*dasin(1.d0)
259  d2r = pi/180.d0
260  re = 6371000.d0
261 
262 c ------------------------------------------------------------------
263 c - Initialize statistical mins and maxes
264 c ------------------------------------------------------------------
265  xlatmin = 90.0000
266  xlatmax = -90.0000
267  xlonmin = 360.0000
268  xlonmax = 0.0000
269 
270 c ------------------------------------------------------------------
271 c - 2016 09 13
272 c - Initialize counts for final report
273 c ------------------------------------------------------------------
274  npts = 0
275  nptslat = 0
276  nptslon = 0
277  nptseht = 0
278 
279 c ------------------------------------------------------------------
280 c - Read the control file in and prepare things.
281 c - For now, presume it is in the exact order below, rather than searching.
282 c ------------------------------------------------------------------
283 c - Header line. Can contain anything.
284  cline='HEADER'
285  read(1,'(a)')card
286  if(card(1:6).ne.cline)then
287  write(6,6000)cline,trim(cfname)
288  stop
289  endif
290  cheader = trim(card(8:200))
291 
292 c - REGION. Must conform to the following list (expand as needed):
293 c conus, alaska, prvi, hawaii, guamcnmi, as, pribilof, stlawrence
294  cline='REGION'
295  read(1,'(a)')card
296  if(card(1:6).ne.cline)then
297  write(6,6000)cline,trim(cfname)
298  stop
299  endif
300  cregion = trim(card(8:200))
301 
302 c - DATUM1. The older datum, chronologically.
303  cline='DATUM1'
304  read(1,'(a)')card
305  if(card(1:6).ne.cline)then
306  write(6,6000)cline,trim(cfname)
307  stop
308  endif
309  cdatum1 = trim(card(8:200))
310 
311 c - DATUM2. The newer datum, chronologically.
312  cline='DATUM2'
313  read(1,'(a)')card
314  if(card(1:6).ne.cline)then
315  write(6,6000)cline,trim(cfname)
316  stop
317  endif
318  cdatum2 = trim(card(8:200))
319 
320 c - REJMET. The rejection criteria in meters.
321 c - Basically if any latitude shift or longitude
322 c - shift or horizontal shift exceeds this value
323 c - (in absolute value), then all shifts for this
324 c - point are set to zero (to avoid asterisks in
325 c - the output file) but the whole line is labeled
326 c - with a triple reject criteria, effectively
327 c - eliminating the pair from use.
328  cline='REJMET'
329  read(1,'(a)')card
330  if(card(1:6).ne.cline)then
331  write(6,6000)cline,trim(cfname)
332  stop
333  endif
334  read(trim(card(8:200)),*)rejmet
335 
336 c - NFILES. The number of *.in files which connect
337 c - the old and new datums in the region being
338 c - addressed.
339  cline='NFILES'
340  read(1,'(a)')card
341  if(card(1:6).ne.cline)then
342  write(6,6000)cline,trim(cfname)
343  stop
344  endif
345  read(card(9:10),*)nfiles
346 
347  6000 format(6x,'makework.f: Expecting ',a6,' line in ',
348  *a,' but not found. Stopping')
349 c ------------------------------------------------------------------
350 c - Read the "manual edits" (workedits) file into RAM, so that it may
351 c - be applied on the fly as we go through the "in" files and build
352 c - the "work" file. Note that the "workedits" file contains every
353 c - manual edit, for every point, for every combination of
354 c - old datum/new datum/region that we might work in. This prevents
355 c - the need for multiple files and allows easy access to every change
356 c - requested.
357 c - 2016 09 13: Update: Also set all "EditTracker" values to ".false."
358 c - so that they can be turned on to ".true." as they actually
359 c - get applied. Then we can look for remainign ".false." ones.
360 c ------------------------------------------------------------------
361 
362  neditstotal = 0
363  neditsrelevant = 0
364  neditsrelevantlat = 0
365  neditsrelevantlon = 0
366  neditsrelevanteht = 0
367 c - 2016 09 13
368  neditsrelevantused = 0
369  neditsrelevantusedlat = 0
370  neditsrelevantusedlon = 0
371  neditsrelevantusedeht = 0
372 
373  701 read(20,702,end=703)card
374 c - If I find a comment or blank, skip out. Otherwise
375 c - pick up one edit's worth of data, and then loop back up
376  if(card(1:1).eq.'#' .or. card(1:1).eq.' ')goto 701
377  neditstotal = neditstotal + 1
378  if(trim(card( 1: 10)) .eq. trim(olddtm) .and.
379  * trim(card( 12: 21)) .eq. trim(newdtm) .and.
380  * trim(card( 23: 32)) .eq. trim(region) )then
381  neditsrelevant = neditsrelevant + 1
382 c - 2016 09 13:
383  edittracker(neditsrelevant) = .false.
384 
385  if(card(41:41).eq.'1')then
386  neditsrelevantlat = neditsrelevantlat + 1
387  endif
388  if(card(42:42).eq.'1')then
389  neditsrelevantlon = neditsrelevantlon + 1
390  endif
391  if(card(43:43).eq.'1')then
392  neditsrelevanteht = neditsrelevanteht + 1
393  endif
394 
395  editolddtm(neditsrelevant) = card( 1: 10)
396  editnewdtm(neditsrelevant) = card( 12: 21)
397  editregion(neditsrelevant) = card( 23: 32)
398  editpid(neditsrelevant) = card( 34: 39)
399  editrejlat(neditsrelevant) = card( 41: 41)
400  editrejlon(neditsrelevant) = card( 42: 42)
401  editrejeht(neditsrelevant) = card( 43: 43)
402 c - 2016 09 13: FORCE the "0" flags to be " "
403  if(editrejlat(neditsrelevant).eq.'0')then
404  editrejlat(neditsrelevant) = ' '
405  endif
406  if(editrejlon(neditsrelevant).eq.'0')then
407  editrejlon(neditsrelevant) = ' '
408  endif
409  if(editrejeht(neditsrelevant).eq.'0')then
410  editrejeht(neditsrelevant) = ' '
411  endif
412  endif
413  goto 701
414  702 format(a)
415 c 703 write(6,704)neditsTotal,neditsRelevant
416 c 704 format(6x,'makework.f: Total Manual Edits Found: ',i6,/,
417 c * 6x,' Relevant Manual Edits Found: ',i6)
418 
419  703 write(6,704)neditstotal,neditsrelevant,
420  *neditsrelevantlat,neditsrelevantlon,neditsrelevanteht
421  704 format(6x,'makework.f: Total Manual Edits Found: ',i6,/,
422  * 6x,' Initial Relevant Manual Edits Found: ',i6,/,
423  * 6x,' ...of these, # in LAT : ',i6,/,
424  * 6x,' # in LON : ',i6,/,
425  * 6x,' # in EHT : ',i6)
426 
427 c - Loop over all *.in files, compute stuff, populate work file.
428 
429  do 1 i=1,nfiles
430  read(1,'(a)')fname0
431  fname='InFiles/'//trim(fname0)
432  write(6,999)trim(fname)
433  999 format(6x,'makework.f: Processing file: ',a)
434  open(10,file=fname,status='old',form='formatted')
435 
436  read(10,100)nameh,namef
437 c write(6,100)nameh,namef
438  2 read(10,101,end=98)pid,state,clath,clonh,cehth,clatf,clonf,cehtf
439  badlat=.false.
440  badlon=.false.
441  badeht=.false.
442  if(clath(11:13).eq.'N/A' .or.
443  * clatf(11:13).eq.'N/A')badlat=.true.
444  if(clonh(12:14).eq.'N/A' .or.
445  * clonf(12:14).eq.'N/A')badlon=.true.
446  if(cehth( 7: 9).eq.'N/A' .or.
447  * cehtf( 7: 9).eq.'N/A')badeht=.true.
448 
449 c - See DRU-10, p. 142 for the below decision
450  if ( badlat .and. badlon .and. badeht)then
451  rejlat = '1'
452  rejlon = '1'
453  rejeht = '1'
454  elseif( badlat .and. badlon .and. .not.badeht)then
455  rejlat = '1'
456  rejlon = '1'
457  rejeht = '2'
458  elseif( badlat .and. .not.badlon .and. badeht)then
459  rejlat = '1'
460  rejlon = '2'
461  rejeht = '1'
462  elseif( badlat .and. .not.badlon .and. .not.badeht)then
463  rejlat = '1'
464  rejlon = '2'
465  rejeht = '2'
466  elseif(.not.badlat .and. badlon .and. badeht)then
467  rejlat = '2'
468  rejlon = '1'
469  rejeht = '1'
470  elseif(.not.badlat .and. badlon .and. .not.badeht)then
471  rejlat = '2'
472  rejlon = '1'
473  rejeht = '2'
474  elseif(.not.badlat .and. .not.badlon .and. badeht)then
475  rejlat = ' '
476  rejlon = ' '
477  rejeht = '1'
478  elseif(.not.badlat .and. .not.badlon .and. .not.badeht)then
479  rejlat = ' '
480  rejlon = ' '
481  rejeht = ' '
482  else
483  stop 10005
484  endif
485 
486 
487  if(.not.badlat)then
488  read(clath(2: 3),'(i2.2)')ilatdh
489  read(clath(4: 5),'(i2.2)')ilatmh
490  read(clath(6:13),'(f8.5)')xlatsh
491  read(clatf(2: 3),'(i2.2)')ilatdf
492  read(clatf(4: 5),'(i2.2)')ilatmf
493  read(clatf(6:13),'(f8.5)')xlatsf
494  xlath = dble(ilatdh) + dble(ilatmh)/60.d0 + xlatsh/3600.d0
495  if(clath(1:1).eq.'S')xlath = -xlath
496  xlatf = dble(ilatdf) + dble(ilatmf)/60.d0 + xlatsf/3600.d0
497  if(clatf(1:1).eq.'S')xlatf = -xlatf
498 
499  dlat = xlatf - xlath
500  dlatsec = dlat * 3600.d0
501  dlatm = dlat*d2r*re
502 
503  endif
504 
505  if(.not.badlon)then
506  read(clonh(2: 4),'(i3.3)')ilondh
507  read(clonh(5: 6),'(i2.2)')ilonmh
508  read(clonh(7:14),'(f8.5)')xlonsh
509  read(clonf(2: 4),'(i3.3)')ilondf
510  read(clonf(5: 6),'(i2.2)')ilonmf
511  read(clonf(7:14),'(f8.5)')xlonsf
512  xlonh = dble(ilondh) + dble(ilonmh)/60.d0 + xlonsh/3600.d0
513  if(clonh(1:1).eq.'W')xlonh = 360.d0 - xlonh
514  xlonf = dble(ilondf) + dble(ilonmf)/60.d0 + xlonsf/3600.d0
515  if(clonf(1:1).eq.'W')xlonf = 360.d0 - xlonf
516 
517 c - Need a cosine of latitude. On the super rare off chance that either the
518 c - HARN or the FBN latitude are "N/A", despite both of the longitudes
519 c - NOT being "N/A", I will look for this eventuality and account for it...
520  if(clath(11:13).ne.'N/A')then
521  coslat = dcos(xlath*d2r)
522  else
523  if(clatf(11:13).ne.'N/A')then
524  coslat = dcos(xlatf*d2r)
525  else
526  coslat = 0.d0
527  rejlon = '3'
528  endif
529  endif
530 
531  dlon = xlonf - xlonh
532  dlonsec = dlon * 3600.d0
533  dlonm = coslat*dlon*d2r*re
534  endif
535 
536  if(.not.badeht)then
537  read(cehth(1: 9),'(f9.3)')xehth
538  read(cehtf(1: 9),'(f9.3)')xehtf
539  dehtm = xehtf - xehth
540  endif
541 
542  if(.not.badlat .and. .not.badlon)then
543  dhorsec = dsqrt(dlatsec**2 + dlonsec**2)
544  dhorm = dsqrt(dlatm**2 + dlonm**2)
545  azhor = datan2(dlonm,dlatm)/d2r
546  if(azhor.lt.0)azhor = azhor + 360.d0
547  endif
548 
549 c - 2015 08 14 -- Do an auto-rejection on ridiculously large shifts
550  if(dabs(dlatm).gt.rejmet .or.
551  * dabs(dlonm).gt.rejmet .or.
552  * dabs(dhorm ).gt.rejmet )then
553  badlat =.true.
554  badlon =.true.
555  badeht =.true.
556  rejlat ='1'
557  rejlon ='1'
558  rejeht ='1'
559  endif
560 
561 c - If there was a bad lat pair, bad lon pair or bad height pair
562 c - then turn this value to 0.0
563  if(badlat)then
564  xlath = 0.d0
565  dlatsec = 0.d0
566  dlatm = 0.d0
567  dhorsec = 0.d0
568  dhorm = 0.d0
569  azhor = 0.d0
570  endif
571 
572  if(badlon)then
573  xlonh = 0.d0
574  dlonsec = 0.d0
575  dlonm = 0.d0
576  dhorsec = 0.d0
577  dhorm = 0.d0
578  azhor = 0.d0
579  endif
580 
581  if(badeht)then
582  xehth = 0.d0
583  dehtm = 0.d0
584  endif
585 
586 
587 c - 2016 09 14: Check if this point is out side of the
588 c - official grid bounds for this region. If so,
589 c - flag it as bad in lat/lon/eht.
590  if(xlath.lt.xs.or.xlath.gt.xn .or.
591  * xlonh.lt.xw.or.xlonh.gt.xe) then
592  noutside = noutside + 1
593  badlat = .true.
594  badlon = .true.
595  badeht = .true.
596  rejlat = '4'
597  rejlon = '4'
598  rejeht = '4'
599  write(6,2003)pid,xlath,xlonh
600  endif
601  2003 format(
602  * 6x,'program makework.f: InFile point found',
603  * ' and flagged with 444 which is outside ',
604  * ' the grid boundaries:',a6,1x,f14.10,1x,f14.10)
605 
606 
607 
608 c - Now, using only GOOD points, let's find the min/max of everything
609  if(.not.badlat .and. .not.badlon)then
610  if(xlath.lt.xlatmin)xlatmin=xlath
611  if(xlath.gt.xlatmax)xlatmax=xlath
612  if(xlonh.lt.xlonmin)xlonmin=xlonh
613  if(xlonh.gt.xlonmax)xlonmax=xlonh
614  endif
615 
616 c - At this point I should be able to print out what I need.
617 c - Cycle through all of the Relevant manual edits, looking
618 c - for this PID, and apply the manual rejections.
619 
620  do 720 irel = 1,neditsrelevant
621 c write(6,751)irel,trim(EditPID(irel)),pid
622  751 format('Checking irel=',i4,' PID=',a6,
623  *' against true pid=',a6)
624  if(trim(editpid(irel)).ne.pid)goto 720
625 
626 c - Updated 2016 09 13: If I get here, I'm using this
627 c - "relevant edit", so track it:
628 c write(6,*) ' HERE'
629  edittracker(irel) = .true.
630 
631 c - Updated 9/13/2016: The final reject code given to "ok" values
632 c - must be a BLANK, not a ZERO! I've switched
633 c - the code earlier, so all "EditRejXXX" values that were zeros
634 c - are blanks. So this code is good here now.
635 c - Modified 2/26/2016: *ONLY* replace a rejection code under these
636 c - rules:
637 c A "0" or " " never replaces a standing "1" (or "2")
638 c A "1" always replaces a standing "0" or " "
639 c (And for ease of coding, a standing "0" or " " can be replaced with another "0" or " ")
640 
641  if(rejlat.eq.'0' .or. rejlat.eq.' ')
642  * rejlat = editrejlat(irel)
643  if(rejlon.eq.'0' .or. rejlon.eq.' ')
644  * rejlon = editrejlon(irel)
645  if(rejeht.eq.'0' .or. rejeht.eq.' ')
646  * rejeht = editrejeht(irel)
647 c rejlat = EditRejLat(irel)
648 c rejlon = EditRejLon(irel)
649 c rejeht = EditRejEht(irel)
650 c - End of 2/26/2016 Modification
651 
652  720 continue
653 
654 c write(6,104)pid,state,rejlat,rejlon,rejeht,xlath,xlonh,xehth,
655 c *dlatsec,dlonsec,dehtm,dhorsec,azhor,dlatm,dlonm,dhorm,
656 c *olddtm,newdtm
657 
658 
659  write(2,104)pid,state,rejlat,rejlon,rejeht,xlath,xlonh,xehth,
660  * dlatsec,dlonsec,dehtm,dhorsec,azhor,dlatm,dlonm,dhorm,
661  * olddtm,newdtm
662 
663 c - 2016 09 13
664  npts = npts + 1
665  if(rejlat.eq.' ')nptslat = nptslat + 1
666  if(rejlon.eq.' ')nptslon = nptslon + 1
667  if(rejeht.eq.' ')nptseht = nptseht + 1
668 
669  104 format(a6,1x,a2,a1,a1,a1,1x,f14.10,1x,f14.10,1x,f8.3,1x,
670  * f9.5,1x,f9.5,1x,f9.3,1x,f9.5,1x,f9.5,1x,f9.3,1x,f9.3,1x,f9.3,
671  * 1x,a10,1x,a10)
672 
673 c - The new file has the following format:
674 c Cols Format Description
675 c 1- 6 a6 PID
676 c 7 1x - blank -
677 c 8- 9 a2 State
678 c 10 a1 Reject code for missing latitude pair (blank for good)
679 c 11 a1 Reject code for missing longitude pair (blank for good)
680 c 12 a1 Reject code for missing ellip ht pair (blank for good)
681 c 13 1x - blank -
682 c 14- 27 f14.10 Latitude (HARN), decimal degrees (-90 to +90)
683 c 28 1x - blank -
684 c 29- 42 f14.10 Lonitude (HARN), decimal degrees (0 to 360)
685 c 43 1x - blank -
686 c 44- 51 f8.3 Ellipsoid Height (HARN), meters
687 c 52 1x - blank -
688 c 53- 61 f9.5 Delta Lat (FBN-HARN), arcseconds
689 c 62 1x - blank -
690 c 63- 71 f9.5 Delta Lon (FBN-HARN), arcseconds
691 c 72 1x - blank -
692 c 73- 81 f9.3 Delta Ell Ht (FBN-HARN), meters
693 c 82 1x - blank -
694 c 83- 91 f9.5 Delta Horizontal (absolute value), arcseconds
695 c 92 1x - blank -
696 c 93-101 f9.5 Azimuth of Delta Horizontal (0-360), degrees
697 c 102 1x - blank -
698 c 103-111 f9.3 Delta Lat (FBN-HARN), meters
699 c 112 1x - blank -
700 c 113-121 f9.3 Delta Lon (FBN-HARN), meters
701 c 122 1x - blank -
702 c 123-131 f9.3 Delta Horizontal (absolute value), meters
703 
704  goto 2
705 
706  98 continue
707 
708 c write(6,103)trim(fname)
709  close(10)
710 c pause
711 
712  1 continue
713 
714 c - Updated 2016 09 13:
715 c - I notice that sometimes there is a typo in the workedits file
716 c - and a "relevant" edit isn't really relevant, because, while
717 c - it matches olddatum/newdatum/region, the actual PID isn't
718 c - even in the InFiles for this olddatum/newdatum/region.
719 c - So in the above code, I've been tracking who got used
720 c - in "EditTracker". Let's spin over all so-called "Relevant Edits"
721 c - and see who did NOT get used.
722 
723 c write(6,*)EditTracker
724 c pause
725 
726  inotrel = 0
727  inotrellat = 0
728  inotrellon = 0
729  inotreleht = 0
730 
731  do 740 irel = 1,neditsrelevant
732  if(edittracker(irel))goto 740
733  inotrel = inotrel + 1
734  if(editrejlat(irel).eq.'1')then
735  inotrellat = inotrellat + 1
736  endif
737  if(editrejlon(irel).eq.'1')then
738  inotrellon = inotrellon + 1
739  endif
740  if(editrejeht(irel).eq.'1')then
741  inotreleht = inotreleht + 1
742  endif
743  write(6,731)editpid(irel)
744  731 format(
745  * 6x,'program makework.f: So-called Relevant Edit not used',
746  * ' since PID is not in the incoming data:',a6)
747  740 continue
748 
749  neditsrelevantused = neditsrelevant - inotrel
750  neditsrelevantusedlat = neditsrelevantlat - inotrellat
751  neditsrelevantusedlon = neditsrelevantlon - inotrellon
752  neditsrelevantusedeht = neditsrelevanteht - inotreleht
753 
754  write(6,732)neditsrelevantused,neditsrelevantusedlat,
755  * neditsrelevantusedlon,neditsrelevantusedeht
756  732 format(6x,'makework.f: ',/,
757  * 6x,' Final Relevant Manual Edits Found: ',i6,/,
758  * 6x,' ...of these, # in LAT : ',i6,/,
759  * 6x,' # in LON : ',i6,/,
760  * 6x,' # in EHT : ',i6)
761 
762 
763 
764  write(6,1005)xlatmin,xlatmax,xlonmin,xlonmax
765  1005 format(
766  *6x,'program makework.f: Minimum latitude: ',f14.10,/,
767  *6x,'program makework.f: Maximum latitude: ',f14.10,/,
768  *6x,'program makework.f: Minimum longitude: ',f14.10,/,
769  *6x,'program makework.f: Maximum longitude: ',f14.10)
770 
771 
772 c - 2016 09 13
773  write(6,1010)npts,nptslat,nptslon,nptseht
774  1010 format(
775  *6x,'program makework.f: Total records in work file: ',i9,/,
776  *6x,'program makework.f: With a usable LAT diff : ',i9,/,
777  *6x,'program makework.f: With a usable LON diff : ',i9,/,
778  *6x,'program makework.f: With a usable EHT diff : ',i9)
779 
780 
781 
782  99 write(6,1003)
783  1003 format('END program makework.f')
784 
785 
786 
787 
788  100 format(27x,a15,26x,a15)
789  101 format(a6,1x,a2,5x,a13,1x,a14,1x,a9,3x,a13,1x,a14,1x,a9)
790  102 format(f15.9,1x,f14.9,1x,f5.1)
791  103 format(6x,'makework.f: Done with file : ',a)
792 
793 
794  end
795  include 'Subs/getgridbounds.f'
makework
program makework
Program to create a work file which will serve as the primary information needed to analyze and creat...
Definition: makework.f:146
getgridbounds
subroutine getgridbounds(region, xn, xs, xw, xe)
Subroutine to collect up the GRID boundaries for use in creating NADCON 5.
Definition: getgridbounds.f:26