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