NADCON5-ng  0.0.1
NADCON5 Next Generation
myrms.f
Go to the documentation of this file.
1 c> \ingroup doers
2 c> Part of the NADCON5 build process, generates `gmtbat05`
3 c>
4 c> Creates a batch file called
5 c>
6 c> gmtbat05.(olddtm).(newdtm).(region).(igridsec)
7 c>
8 c>
9 c> Program to
10 c> 1. Run a customized RMS-computing
11 c> algorithm on vector differences (gridded minus true)
12 c> 2. Save the RMS data to GMT-ready plottable files
13 c> 3. Create a GMT batch file to plot both the
14 c> thinned data and removed data in both coverage and vectors
15 c>
16 c> Unlike "mymedian5.f", this program is
17 c> set up to compute the RMS of vector differences
18 c> in a cell-by-cell basis (aka NOT a median filter
19 c> at all, but a true RMS representation of of disagreements)
20 c>
21 c> A "value" is the differential vector of:
22 c>
23 c> interpolated-from-grid minus true
24 c>
25 c> For any cell with at least ONE point with a value, the following is done:
26 c> 1. Compute the average latitude of all points in the cell
27 c> 2. Compute the average longitude of all points in the cell
28 c> 3. Compute the RMS of all values in the cell
29 c>
30 c> The output vector will then reflect these three values.
31 c>
32 c> For latitude and ellipsoid height, the azimuth of the vector
33 c> will ALWAYS be 0.0 (pointing up) while for longitude it will
34 c> always be 90.0 (pointing right). However, these are mere conventions
35 c> as they are not directional vectors anyway, but rather quanta which
36 c> will be gridded and it is the grid which is of utmost importance.
37 c>
38 c> No PIDS will be in the output files, as the output RMS vectors are
39 c> not reflective of any one point, but rather a cell-wide conglomeration
40 c> of information.
41 c>
42 c> See DRU-11, p. 130
43 c>
44 c> ### Program arguments
45 c> Arguments are newline terminated and read from standard input
46 c>
47 c> They are enumerated here
48 c> \param oldtm Source Datum
49 c> \param newdtm Target Datum,region
50 c> \param region Conversion Region
51 c> \param agridsec Grid Spacing in arcsec
52 c>
53 c> Example:
54 c>
55 c> olddatum = 'ussd'
56 c> newdatum = 'nad27'
57 c> region = 'conus'
58 c> agridsec = '900'
59 c>
60 c> ### Program Inputs:
61 c>
62 c> ## Changelog
63 c>
64 c> ### 2016 01 21:
65 c> Updated to RETURN to an old way of registering RMS vectors
66 c> at AVE lat/lon rather than center of cell.
67 c>
68 c> ### 2015 10 28
69 c> Updated to work with new naming scheme (see DRU-11, p. 150) and
70 c> to adopt the central lat/lon for the RMS vectors, rather than
71 c> ave lat/ave lon (see DRU-11, p.145), and to also set up the gridding
72 c> of RMS vectors at the T=0.9 level (see DRU-11, p. 148)
73 c>
74 c> Also added "donzd" functionality to help control the magnitude
75 c> of the Length of Reference Vector on Ground variables.
76 c>
77 c> Also, added a section at the end to create the TOTAL error grids.
78 c> by RMS-combining the "method noise grid" (the "d3" grid)
79 c> with the "data noise grid" (the "rdd...09" grid) into
80 c> one single "transformation error grid"
81 c>
82 c> ### 2016 08 25:
83 c> For reasons that are difficult to describe, "donzd" is
84 c> now "onzd2.f" in /home/dru/Subs. Change and recompile...
85 c>
86 c> ### 2015 10 09:
87 c> Updated to add HOR vectors
88 c>
89 c> ### 2015 10 06:
90 c> Updated
91 c>
92 c> ### 2015 09 16:
93 c> Initial Release, For use in creating NADCON5
94 c> Built by Dru Smith
95 c>
96 c>
97  program myrms
98 
99  implicit real*8(a-h,o-z)
100 
101 c - Allow up to 280,000 points to come in
102 c - which also translates to letting 280,000 grid cells
103 c - to be non-empty.
104  parameter(max=280000)
105 
106 c - And using a little hit or miss logic, the number
107 c - of possible points in one cell can't be more than
108 c - about 5000, right?
109  parameter(maxppcell=5000)
110 
111 c - Holding place for each record, in read-order
112  character*80 cards(max)
113  real*8 xlat(max),xlon(max),z0(max),zm(max)
114 
115 c - Index of what latititude row and longitude column
116 c - is associated with each non-empty cell, in
117 c - cell-found order.
118 c - 1 <= ilapcell() <= nla
119 c - 1 <= ilopcell() <= nlo
120  integer*4 ilapcell(max),ilopcell(max)
121  integer*4 poscode(max)
122  integer*4 indx(max)
123 c - Distance of a point from its cell's center
124  real*8 dist(max)
125 
126 c - Index telling me which cell a given point is
127 c - associated with, in read-order
128 c - 1 <= whichcell() <= ncells
129  integer*4 whichcell(max)
130 
131 c - ppcell(i) contains the count of how many values fall
132 c - in cell "i", where "i" ratchets up from 1 to max
133 c - every time I find a new non-empty cell.
134  integer*4 ppcell(max)
135 
136  character*10 olddtm,newdtm,region
137  character*200 suffix1,suffix2
138  character*200 suffix2t04,suffix2t09,suffix2d3
139 
140 c - Input GMT-ready differentical vector files containing all points:
141  character*200 gfnvmaddlat,gfnvmaddlon,gfnvmaddeht
142  character*200 gfnvsaddlat,gfnvsaddlon
143  character*200 gfnvsaddhor,gfnvmaddhor ! 2015 10 09
144 c - Output GMT-ready RMS differential vector files:
145  character*200 gfnvmrddlat,gfnvmrddlon,gfnvmrddeht
146  character*200 gfnvsrddlat,gfnvsrddlon
147  character*200 gfnvsrddhor,gfnvmrddhor ! 2015 10 09
148 c - Output RMS GMT(surface)-ready RMS differential vector files:
149  character*200 gfnsmrddlat,gfnsmrddlon,gfnsmrddeht
150  character*200 gfnssrddlat,gfnssrddlon
151 c - Output RMS GMT-read RMS differential vector COVERAGE files:
152  character*200 gfncvrddlat,gfncvrddlon,gfncvrddeht
153 c - To be created ".grd" files by surface
154  character*200 rfnvmrddlat,rfnvmrddlon,rfnvmrddeht
155  character*200 rfnvsrddlat,rfnvsrddlon
156 c - To be created ".xyz" files by GMT's grd2xyz program
157  character*200 zfnvmrddlat,zfnvmrddlon,zfnvmrddeht
158  character*200 zfnvsrddlat,zfnvsrddlon
159 c - To be created ".b" files by grd2b
160  character*200 bfnvmrddlat,bfnvmrddlon,bfnvmrddeht
161  character*200 bfnvsrddlat,bfnvsrddlon
162 c - Stats file name
163  character*200 sfn
164  logical*1 nothinlat,nothinlon,nothineht
165 
166 c - "Method Noise" grids:
167  character*200 sadbfnvmtcdlat1000,sadbfnvmtcdlon1000
168  character*200 sadbfnvmtcdeht1000
169  character*200 sadbfnvstcdlat1000,sadbfnvstcdlon1000
170 c -" Total Error" grids in ".b" format:
171  character*200 bfnvsetelat,bfnvsetelon
172  character*200 bfnvmetelat,bfnvmetelon
173  character*200 bfnvmeteeht
174 c -" Total Error" grids in "grd" format:
175  character*200 rfnvsetelat,rfnvsetelon
176  character*200 rfnvmetelat,rfnvmetelon
177  character*200 rfnvmeteeht
178 
179 
180 c - Dummy
181  character*200 fused,gfndv,gfnsu
182 c - Output GMT-batch file:
183  character*200 gmtfile
184 c - Grid spacing, converted to characters for file names
185 c - Presumes:
186 c - 1) igridsec is an integer
187 c - 2) 1 <= igridsec <= 99999
188  integer*4 igridsec
189  character*5 agridsec
190 
191  character*6 pid
192  character*80 card
193 
194  real*8 lorvopc,lorvog0,lorvogm
195  real*8 lorvoghorddm,lorvoghordds
196 
197  real*8 basedd(5)
198 
199 c --------------------------------------
200 c --------------------------------------
201 c --------------------------------------
202 c -------BEGIN GUTS OF PROGRAM----------
203 c --------------------------------------
204 c --------------------------------------
205 c --------------------------------------
206  write(6,1001)
207  1001 format('BEGIN program myrms.f')
208 
209 c ------------------------------------------------------------------
210 c - Some required constants.
211 c ------------------------------------------------------------------
212  lorvopc = 1.d0
213  pi = 2.d0*dasin(1.d0)
214  d2r = pi/180.d0
215  re = 6371000.d0
216  s2m = (1/3600.d0)*d2r*re
217  multiplierlorvog = 2
218 
219 c ------------------------------------------------------------------
220 c - Read in arguments from batch file
221 c ------------------------------------------------------------------
222  read(5,'(a)')olddtm
223  read(5,'(a)')newdtm
224  read(5,'(a)')region
225  read(5,'(a)')agridsec
226 
227 c ------------------------------------------------------------------
228 c - Generate the suffixes used in all our files
229 c ------------------------------------------------------------------
230  read(agridsec,*)igridsec
231  suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region)
232  suffix2=trim(suffix1)//'.'//trim(agridsec)
233  suffix2t04=trim(suffix2)//'.04'
234  suffix2t09=trim(suffix2)//'.09'
235  suffix2d3=trim(suffix2)//'.d3'
236 
237 c ------------------------------------------------------------------
238 c - Open the GMT batch file for gridding the RMS data
239 c ------------------------------------------------------------------
240  gmtfile = 'gmtbat05.'//trim(suffix2)
241  open(99,file=gmtfile,status='new',form='formatted')
242  write(6,1011)trim(gmtfile)
243  1011 format(6x,'myrms.f: Creating GMT batch file ',a)
244  write(99,1030)trim(gmtfile)
245  1030 format('echo BEGIN batch file ',a)
246 
247 
248 c ------------------------------------------------------------------
249 c - Open the "dvstats" file generated by "checkgrid", and use it
250 c - to determine which, if any, elements (lat, lon, eht) to skip
251 c ------------------------------------------------------------------
252  sfn = 'dvstats.'//trim(suffix2)
253  open(2,file=sfn,status='old',form='formatted')
254  write(6,1003)trim(sfn)
255  1003 format(6x,'myrms.f: Opening exisiting stats file ',a)
256  read(2,*)nlat,rmslat
257  read(2,*)nlon,rmslon
258  read(2,*)neht,rmseht
259  nothinlat = .false.
260  nothinlon = .false.
261  nothineht = .false.
262  if(nlat.eq.0)nothinlat=.true.
263  if(nlon.eq.0)nothinlon=.true.
264  if(neht.eq.0)nothineht=.true.
265 
266 c ------------------------------------------------------------------
267 c - Open the input, raw, differential vector GMT-ready files
268 c - Note, unlike "mymedian5", these file names already
269 c - contain the grid spacing, since they refer to differences
270 c - relative to a grid that is grid-spacing dependent.
271 c ------------------------------------------------------------------
272  if(.not.nothinlat)then
273  gfnvsaddlat = 'vsaddlat.'//trim(suffix2)
274  open(21,file=gfnvsaddlat,status='old',form='formatted')
275  write(6,1010)trim(gfnvsaddlat)
276 
277  gfnvmaddlat = 'vmaddlat.'//trim(suffix2)
278  open(26,file=gfnvmaddlat,status='old',form='formatted')
279  write(6,1010)trim(gfnvmaddlat)
280  endif
281 
282  if(.not.nothinlon)then
283  gfnvsaddlon = 'vsaddlon.'//trim(suffix2)
284  open(22,file=gfnvsaddlon,status='old',form='formatted')
285  write(6,1010)trim(gfnvsaddlon)
286 
287  gfnvmaddlon = 'vmaddlon.'//trim(suffix2)
288  open(27,file=gfnvmaddlon,status='old',form='formatted')
289  write(6,1010)trim(gfnvmaddlon)
290  endif
291 
292  if(.not.nothineht)then
293  gfnvmaddeht = 'vmaddeht.'//trim(suffix2)
294  open(23,file=gfnvmaddeht,status='old',form='formatted')
295  write(6,1010)trim(gfnvmaddeht)
296  endif
297 
298  if(.not.nothinlat .and. .not.nothinlon)then
299  gfnvsaddhor = 'vsaddhor.'//trim(suffix2)
300  open(24,file=gfnvsaddhor,status='old',form='formatted')
301  write(6,1010)trim(gfnvsaddhor)
302 
303  gfnvmaddhor = 'vmaddhor.'//trim(suffix2)
304  open(25,file=gfnvmaddhor,status='old',form='formatted')
305  write(6,1010)trim(gfnvmaddhor)
306  endif
307 
308  1010 format(6x,'myrms.f: ',
309  *'Opening existing raw differential vector file ',a)
310 
311 c ------------------------------------------------------------------
312 c - Open the output, RMS'd, GMT-ready DIFFERENTIAL VECTOR files
313 c - "r" means "rms of differential"
314 c ------------------------------------------------------------------
315  if(.not.nothinlat)then
316  gfnvsrddlat = 'vsrddlat.'//trim(suffix2)
317  open(41,file=gfnvsrddlat,status='new',form='formatted')
318  write(6,1012)trim(gfnvsrddlat)
319 
320  gfnvmrddlat = 'vmrddlat.'//trim(suffix2)
321  open(46,file=gfnvmrddlat,status='new',form='formatted')
322  write(6,1012)trim(gfnvmrddlat)
323  endif
324 
325  if(.not.nothinlon)then
326  gfnvsrddlon = 'vsrddlon.'//trim(suffix2)
327  open(42,file=gfnvsrddlon,status='new',form='formatted')
328  write(6,1012)trim(gfnvsrddlon)
329 
330  gfnvmrddlon = 'vmrddlon.'//trim(suffix2)
331  open(47,file=gfnvmrddlon,status='new',form='formatted')
332  write(6,1012)trim(gfnvmrddlon)
333  endif
334 
335  if(.not.nothineht)then
336  gfnvmrddeht = 'vmrddeht.'//trim(suffix2)
337  open(43,file=gfnvmrddeht,status='new',form='formatted')
338  write(6,1012)trim(gfnvmrddeht)
339  endif
340 
341  if(.not.nothinlat .and. .not.nothinlon)then
342  gfnvsrddhor = 'vsrddhor.'//trim(suffix2)
343  open(44,file=gfnvsrddhor,status='new',form='formatted')
344  write(6,1010)trim(gfnvsrddhor)
345 
346  gfnvmrddhor = 'vmrddhor.'//trim(suffix2)
347  open(45,file=gfnvmrddhor,status='new',form='formatted')
348  write(6,1010)trim(gfnvmrddhor)
349  endif
350 
351  1012 format(6x,'myrms.f: ',
352  *'Opening output RMS differential vector file ',a)
353 
354 c ------------------------------------------------------------------
355 c - Open the output, RMS'd, GMT-ready differential COVERAGE files
356 c - "r" means "rms of differential"
357 c ------------------------------------------------------------------
358  if(.not.nothinlat)then
359  gfncvrddlat = 'cvrddlat.'//trim(suffix2)
360  open(51,file=gfncvrddlat,status='new',form='formatted')
361  write(6,1013)trim(gfncvrddlat)
362  endif
363 
364  if(.not.nothinlon)then
365  gfncvrddlon = 'cvrddlon.'//trim(suffix2)
366  open(52,file=gfncvrddlon,status='new',form='formatted')
367  write(6,1013)trim(gfncvrddlon)
368  endif
369 
370  if(.not.nothineht)then
371  gfncvrddeht = 'cvrddeht.'//trim(suffix2)
372  open(53,file=gfncvrddeht,status='new',form='formatted')
373  write(6,1013)trim(gfncvrddeht)
374  endif
375 
376  1013 format(6x,'myrms.f: ',
377  *'Opening output RMS differential coverage file ',a)
378 
379 c ------------------------------------------------------------------
380 c - Open the output, RMS'd GMT-ready files for pushing
381 c - through "surface" to get a grid.
382 c - "r" meaning "RMS'd"
383 c ------------------------------------------------------------------
384  if(.not.nothinlat)then
385  gfnssrddlat = 'ssrddlat.'//trim(suffix2)
386  open(71,file=gfnssrddlat,status='new',form='formatted')
387  write(6,1014)trim(gfnssrddlat)
388 
389  gfnsmrddlat = 'smrddlat.'//trim(suffix2)
390  open(76,file=gfnsmrddlat,status='new',form='formatted')
391  write(6,1014)trim(gfnsmrddlat)
392  endif
393 
394  if(.not.nothinlon)then
395  gfnssrddlon = 'ssrddlon.'//trim(suffix2)
396  open(72,file=gfnssrddlon,status='new',form='formatted')
397  write(6,1014)trim(gfnssrddlon)
398 
399  gfnsmrddlon = 'smrddlon.'//trim(suffix2)
400  open(77,file=gfnsmrddlon,status='new',form='formatted')
401  write(6,1014)trim(gfnsmrddlon)
402  endif
403 
404  if(.not.nothineht)then
405  gfnsmrddeht = 'smrddeht.'//trim(suffix2)
406  open(73,file=gfnsmrddeht,status='new',form='formatted')
407  write(6,1014)trim(gfnsmrddeht)
408  endif
409 
410  1014 format(6x,'myrms.f: ',
411  *'Opening output RMS differential vector ',
412  *'file (for surface):',a)
413 
414 c ------------------------------------------------------------------
415 c - Each region has officially chosen boundaries (which may
416 c - or may not agree with the MAP boundaries). Get the
417 c - official grid boundaries here. See DRU-11, p. 126
418 c ------------------------------------------------------------------
419  call getgridbounds(region,glamx,glamn,glomn,glomx)
420 
421  write(6,1004)trim(region),glamn,glamx,glomn,glomx
422  1004 format(6x,'myrms.f: Region= ',a,/,
423  * 6x,'myrms.f: North = ',f12.6,/,
424  * 6x,'myrms.f: South = ',f12.6,/,
425  * 6x,'myrms.f: West = ',f12.6,/,
426  * 6x,'myrms.f: East = ',f12.6)
427 
428 c -------------------------------------------------------
429 c - Get the header information necessary for a ".b" file
430 c - 1) Convert "igridsec" to decimal degrees
431 c - 2) Count rows and columns
432 c -------------------------------------------------------
433  dgla = dble(igridsec)/3600.d0
434  dglo = dble(igridsec)/3600.d0
435 
436  nla=idnint((glamx-glamn)/dgla)+1
437  nlo=idnint((glomx-glomn)/dglo)+1
438 
439  write(6,3001)glamn,glamx,glomn,glomx,dgla,dglo,nla,nlo
440  3001 format(6x,'myrms.f Cell Structure:',/,
441  *8x,'North = ',f16.10,/,
442  *8x,'South = ',f16.10,/,
443  *8x,'West = ',f16.10,/,
444  *8x,'East = ',f16.10,/,
445  *8x,'DLat = ',f16.10,/,
446  *8x,'DLon = ',f16.10,/,
447  *8x,'NLat = ',i16 ,/,
448  *8x,'NLon = ',i16 )
449 
450 c glamx=glamn+dgla*(nla-1) !*** register the far boundary
451 c glomx=glomn+dglo*(nlo-1)
452 
453 c -----------------------------------------------------------------------
454 c -----------------------------------------------------------------------
455 c -----------------------------------------------------------------------
456 c -----------------------------------------------------------------------
457 c - MAIN LOOP
458 c - First trip is to process on differential
459 c - latitude vectors (both units). Second trip
460 c - on differential longitude vectors (both units).
461 c - Third trip on differential ellipsoid height vectors (meters only).
462 c - What's done in each trip:
463 c - 1) Data are sorted into respective cells
464 c - 2) An overall RMS value for the whole dataset is tallied
465 c - 3) After all data are in RAM, each cell's average latitude,
466 c - average longitude and RMS value are computed
467 c - 4) The data in #3 is put (including proper scaling, based on
468 c - the overall RMS value) into a vector plot file
469 c - 5) The data in #3 is put (no scaling necessary) into a
470 c - file ready to go into surface for gridding.
471 c -----------------------------------------------------------------------
472 c -----------------------------------------------------------------------
473 c -----------------------------------------------------------------------
474 c -----------------------------------------------------------------------
475 
476  nrmslat = 0
477  nrmslon = 0
478  nrmseht = 0
479 
480  do 2001 iloop = 1,3
481 
482  if(iloop.eq.1)then
483  write(6,2002)trim(gfnvsrddlat)
484  write(6,2002)trim(gfnvmrddlat)
485  endif
486  if(iloop.eq.2)then
487  write(6,2002)trim(gfnvsrddlon)
488  write(6,2002)trim(gfnvmrddlon)
489  endif
490  if(iloop.eq.3)then
491  write(6,2002)trim(gfnvmrddeht)
492  endif
493 
494  if(iloop.eq.1 .and. nothinlat)goto 2001
495  if(iloop.eq.2 .and. nothinlon)goto 2001
496  if(iloop.eq.3 .and. nothineht)goto 2001
497 
498  lin = 20 + iloop
499 
500  if(iloop.eq.1 .or. iloop.eq.3)az=0.d0
501  if(iloop.eq.2 )az=90.d0
502 
503  2002 format(6x,'myrms.f : RMS computing on file :',a)
504 
505 
506 c 101 format(f16.10,1x,f15.10,1x,6x ,1x,12x ,1x,5x ,1x,f9.5,1x,a6)
507 c - To pull arcseconds:
508  101 format(f16.10,1x,f15.10,1x, 6x,1x, 12x,1x,f9.5,1x, 9x,1x,a6)
509 c - To pull meters:
510  1101 format(f16.10,1x,f15.10,1x, 6x,1x, 12x,1x, 9x,1x,f9.3,1x,a6)
511 
512 
513 
514  2003 format(6x,'myrms.f : Point outside boundaries: ',
515  * f16.10,1x,f15.10,1x,f9.5,1x,a6)
516 
517 c - Set the count of how many values are in each non-empty
518 c - grid cell to zero, as well as the count of how many cells
519 c - have data to zero
520  do 201 i=1,max
521  ppcell(i)=0
522  201 continue
523  ncells = 0
524 
525 c - Top of READ loop
526 c - "ikt" is the count of all records in the vector file
527 c - whether we keep them or not. "ipid" is the count
528 c - points found that are inside our grid boundaries.
529  ikt = 0
530  ipid = 0
531 
532 c - rmsoverallm is the RMS of all vectors in one file,
533 c - converted into meters.
534  rmsoverallm = 0.d0
535 
536  100 read(lin,'(a)',end=777)card
537  if(iloop.le.2)read(card, 101)glo,gla,z,pid
538  if(iloop.eq.3)read(card,1101)glo,gla,z,pid
539  ikt = ikt + 1
540 
541 c - Store the read data and values both
542  cards(ikt) = card
543 c - z0 = primary coordinates
544  z0(ikt) = z
545  xlat(ikt) = gla
546  xlon(ikt) = glo
547  if (iloop.eq.1)then
548  zm(ikt) = z * s2m
549  elseif(iloop.eq.2)then
550  coslat = dcos(gla*d2r)
551  zm(ikt) = z * s2m * coslat
552  else
553  zm(ikt)= z
554  endif
555  rms0 = rms0 + z0(ikt)**2
556  rmsm = rmsm + zm(ikt)**2
557 
558 c - Double check that the data we're using is inside our
559 c - pre-arranged grid boundaries for this region
560  if(gla.lt.glamn.or.gla.gt.glamx .or.
561  * glo.lt.glomn.or.glo.gt.glomx) then
562  write(6,2003)gla,glo,z,pid
563  goto 100
564  endif
565 
566  ipid = ipid + 1
567 
568 c - Determine which row and column our point falls in. Then
569 c - convert that combination of row/col into a single
570 c - value (encoded as ila*100000+ilo)
571  ila=idnint((gla-glamn)/dgla)+1
572  ilo=idnint((glo-glomn)/dglo)+1
573  ipos=ila*100000+ilo
574 
575 c - Is this point a new cell?
576  do 202 j=1,ncells
577 
578 c - OLD cell
579 c - Ratchet up the count of points in this cell
580 c - Store the unique input number as being in this cell
581 c - Jump out
582  if(ipos.eq.poscode(j))then
583  ppcell(j) = ppcell(j) + 1
584  whichcell(ikt) = j
585  goto 203
586  endif
587  202 continue
588 
589 c - NEW cell (jump over this code from previous loop if
590 c - an old cell were found)
591 
592  ncells = ncells + 1
593  poscode(ncells) = ipos
594  ilapcell(ncells) = ila
595  ilopcell(ncells) = ilo
596  ppcell(ncells) = 1
597  whichcell(ikt) = ncells
598 
599  203 continue
600 
601  goto 100
602 
603  777 continue
604 
605 c - At this point, I've read ONE differential vector
606 c - file into RAM and identified the cell of each
607 c - entry.
608 
609  nkt = ikt
610  npid = ipid
611 
612  rms0 = dsqrt(rms0/nkt)
613  rmsm = dsqrt(rmsm/nkt)
614 
615  if(iloop.eq.1)nrmslat = nkt
616  if(iloop.eq.2)nrmslon = nkt
617  if(iloop.eq.3)nrmseht = nkt
618 
619 c - Use the RMS value to compute scales needed for our output vector file.
620 c - Why not AVE? Because the AVE is expected to be VERY close to zero
621 c - but RMS represents a better overall level of disagreement.
622 
623 c iq0 = floor(log10(rms0))
624 c qq0 = 10.d0**(iq0-1)
625 c lorvog0 = MultiplierLorvog*qq0*nint(rms0/qq0)
626  lorvog0 = onzd2(multiplierlorvog*rms0)
627  g02pc = lorvopc / lorvog0
628 
629 c iqm = floor(log10(rmsm))
630 c qqm = 10.d0**(iqm-1)
631 c lorvogm = MultiplierLorvog*qqm*nint(rmsm/qqm)
632  lorvogm = onzd2(multiplierlorvog*rmsm)
633  gm2pc = lorvopc / lorvogm
634 
635 c - 2015 10 09
636  if(iloop.eq.1)then
637  basedd(1) = rms0
638  basedd(4) = rmsm
639  endif
640  if(iloop.eq.2)then
641  basedd(2) = rms0
642  basedd(5) = rmsm
643  endif
644 
645 
646  write(6,2004) nkt,npid,igridsec,ncells
647  2004 format(
648  *6x,'myrms.f: Done with file ',/,
649  *6x,'myrms.f: Points in File : ',i10,/,
650  *6x,'myrms.f: Points within Grid Bounds : ',i10,/,
651  *6x,'myrms.f: Cell Size (Arcseconds) : ',i10,/,
652  *6x,'myrms.f: Number of Cells with Data : ',i10)
653  write(6,2005)
654  2005 format(6x,'myrms.f: Begin RMS computations')
655 
656 c - Sort all of our data (nkt values) by their poscode
657  write(6,2020)
658  2020 format(6x,'myrms.f : Sorting data...')
659  call indexxi(nkt,max,whichcell,indx)
660 
661 c - Spit out our data in WHICHCELL order
662 c do 400 i=1,nkt
663 c j = indx(i)
664 c icell = whichcell(j)
665 c if(mod(i,1000).eq.0)then
666 c write(6,2010)cards(j),icell,ilapcell(icell),
667 c * ilopcell(icell),poscode(icell)
668 c endif
669 c 400 continue
670  2010 format(i8,1x,a80,1x,4(i10))
671 
672 c - 2015/11/06: Change from average lat/lon to
673 c - center of cell for registration of RMS vector,
674 c - as per discussions with Andria
675 
676 c - Now go through each cell and compute
677 c - the average latitude, average longitude
678 c - and RMS value. As each alat/alon/rval is found,
679 c - put it in the appropriate output vector file.
680 c - Change as of 10/28/2015: No longer using
681 c - average lat/lon, but rather center of cell.
682 
683  inumlo = 0
684  inumhi = 0
685  ikt = 0
686 
687  if(iloop.eq.1)then
688  write(6,2021)trim(gfnvsrddlat),trim(gfnssrddlat)
689  write(6,2021)trim(gfnvmrddlat),trim(gfnsmrddlat)
690  endif
691  if(iloop.eq.2)then
692  write(6,2021)trim(gfnvsrddlon),trim(gfnssrddlon)
693  write(6,2021)trim(gfnvmrddlon),trim(gfnsmrddlon)
694  endif
695  if(iloop.eq.3)then
696  write(6,2021)trim(gfnvmrddeht),trim(gfnsmrddeht)
697  endif
698 
699  2021 format(6x,'myrms.f : ',
700  * 'Populating RMS diff vec file: ',a,/,
701  * 6x,'myrms.f : ',
702  * 'Populating RMS diff vec surface-ready file: ',a)
703 
704  ifileout1 = 40+iloop ! DD vectors, primary units
705  ifileout3 = 50+iloop ! DD coverage
706  ifileout2 = 70+iloop ! DD vectors for surface, primary units
707 
708  do 401 icell=1,ncells
709  inumlo = inumhi + 1
710  inumhi = inumlo + ppcell(icell) - 1
711 
712 c - 10/28/2015: removed...useless, right?
713 c do 402 jpt=1,ppcell(icell)
714 c ikt = ikt + 1
715 c j = indx(ikt)
716 c 402 continue
717 
718  call getrms(z0,max,inumlo,inumhi,indx,
719  * maxppcell,xlat,xlon,alat,alon,rval0)
720 
721 c - 01/21/2016: Returned to use of AVE lat/lon for registering RMS vectors
722 c - 10/28/2015: Use center of cell
723 c ila = ilapcell(icell)
724 c ilo = ilopcell(icell)
725 c glacc = glamn + (ila-1)*dgla + (dgla/2.d0)
726 c glocc = glomn + (ilo-1)*dglo + (dglo/2.d0)
727 c alat = glacc
728 c alon = glocc
729 
730  if (iloop.eq.1)then
731  rvalm = rval0 * s2m
732  elseif(iloop.eq.2)then
733  coslat = dcos(alat*d2r)
734  rvalm = rval0 * s2m * coslat
735  else
736  rvalm = rval0
737  endif
738  vc0 = dabs(rval0 * g02pc)
739  vcm = dabs(rvalm * gm2pc)
740 
741 c - For vector plotting:
742  write(ifileout1 ,2101)alon,alat,az,vc0,rval0,rvalm
743  if(iloop.lt.3)
744  * write(ifileout1+5,2101)alon,alat,az,vcm,rval0,rvalm
745  2101 format(f16.10,1x,f15.10,1x,f6.2,1x,
746  * f12.2,1x,f9.5,1x,f9.3)
747 
748 c - For coverage plotting:
749  write(ifileout3,2103)alon,alat,1.0
750  2103 format(f16.10,1x,f15.10,1x,f6.2)
751 
752 c - For use in surface for gridding:
753  write(ifileout2 ,2102)alon,alat,rval0
754  if(iloop.lt.3)
755  * write(ifileout2+5,2102)alon,alat,rvalm
756  2102 format(f16.10,1x,f15.10,1x,f10.5)
757 
758  401 continue
759 
760  2001 continue
761 
762 c -------------------------------------------------
763 c - 2015 10 09
764 c - Spin through the newly created RMS DD vectors in
765 c - lat/lon (in sec and meters) and create an RMS DD
766 c - vector file in HOR (in sec and meters)
767 c
768 c - Do NOT create this from the "add" horizontal
769 c - files, but rather from the newly created "rdd"
770 c - lat and lon files.
771 c -------------------------------------------------
772 
773  baseddhors = sqrt(basedd(1)**2+basedd(2)**2)
774  baseddhorm = sqrt(basedd(4)**2+basedd(5)**2)
775 
776 c iqhorddm = floor(log10(baseddhorm))
777 c qqhorddm = 10.d0**(iqhorddm-1)
778 c lorvoghorddm = MultiplierLorvog*qqhorddm*nint(baseddhorm/qqhorddm)
779  lorvoghorddm = onzd2(multiplierlorvog*baseddhorm)
780  gm2pchordd = lorvopc / lorvoghorddm
781 
782 c iqhordds = floor(log10(baseddhors))
783 c qqhordds = 10.d0**(iqhordds-1)
784 c lorvoghordds = MultiplierLorvog*qqhordds*nint(baseddhors/qqhordds)
785  lorvoghordds = onzd2(multiplierlorvog*baseddhors)
786  gs2pchordd = lorvopc / lorvoghordds
787 
788  write(6,2504)
789  2504 format
790  *(6x,'myrms.f: Populating RMS dd horizontal vector files')
791  iflats = 41
792  iflons = 42
793  iflatm = 46
794  iflonm = 47
795  ifhors = 44
796  ifhorm = 45
797  rewind(iflats)
798  rewind(iflons)
799  rewind(iflatm)
800  rewind(iflonm)
801  2503 read(iflats,2101,end=2502)xlo1,xla1,az1,vc1,xs1,xm1
802  read(iflons,2101)xlo2,xla2,az2,vc2,xs2,xm2
803  if(xlo1.ne.xlo2)stop 10501
804  if(xla1.ne.xla2)stop 10502
805 c if(az1.ne.90.00)stop 10503
806 c if(az2.ne. 0.00)stop 10504
807  azhor = datan2(xm2,xm1)/d2r
808  if(azhor.lt.0)azhor = azhor + 360.d0
809  xshor = dsqrt(xs1**2 + xs2**2)
810  xmhor = dsqrt(xm1**2 + xm2**2)
811  vchors = xshor * gs2pchordd
812  vchorm = xmhor * gm2pchordd
813  write(ifhors,2101)xlo1,xla1,azhor,vchors,xshor,xmhor
814  write(ifhorm,2101)xlo1,xla1,azhor,vchorm,xshor,xmhor
815  goto 2503
816  2502 continue
817 
818 c -------------------------------------------------------
819 c - Based on everything we have so far, put a bunch
820 c - of "surface" calls into the GMT batch file.
821 c - These will be to turn the RMS'd differential vector data
822 c - into grids. DO NOT GENERATE A CALL TO SURFACE
823 c - IF THE NUMBER OF RMS'd DATA IS ZERO.
824 
825 c - Because "surface", as part of GMT, returns its
826 c - own binary grid format (called ".grd" herein), it must be
827 c - converted to ".b" format. The easiest way to
828 c - do so is to use the GMT built-in routine "grd2xyz"
829 c - with the "-bos" extension to create a binary XYZ
830 c - list file. Then run Dennis's homemade "xyz2b.for"
831 c - routine to finally arrive at the ".b" binary grid format.
832 c
833 c - 2015/11/06: Updated to put the ".09." in the names
834 c -------------------------------------------------------
835  cmidlat=dcos(((glamn+glamx)/2.d0)*d2r)
836  rfnvsrddlat = 'vsrddlat.'//trim(suffix2t09)//'.grd'
837  rfnvsrddlon = 'vsrddlon.'//trim(suffix2t09)//'.grd'
838  rfnvmrddeht = 'vmrddeht.'//trim(suffix2t09)//'.grd'
839  rfnvmrddlat = 'vmrddlat.'//trim(suffix2t09)//'.grd'
840  rfnvmrddlon = 'vmrddlon.'//trim(suffix2t09)//'.grd'
841 
842  zfnvsrddlat = 'vsrddlat.'//trim(suffix2t09)//'.xyz'
843  zfnvsrddlon = 'vsrddlon.'//trim(suffix2t09)//'.xyz'
844  zfnvmrddeht = 'vmrddeht.'//trim(suffix2t09)//'.xyz'
845  zfnvmrddlat = 'vmrddlat.'//trim(suffix2t09)//'.xyz'
846  zfnvmrddlon = 'vmrddlon.'//trim(suffix2t09)//'.xyz'
847 
848  bfnvsrddlat = 'vsrddlat.'//trim(suffix2t09)//'.b'
849  bfnvsrddlon = 'vsrddlon.'//trim(suffix2t09)//'.b'
850  bfnvmrddeht = 'vmrddeht.'//trim(suffix2t09)//'.b'
851  bfnvmrddlat = 'vmrddlat.'//trim(suffix2t09)//'.b'
852  bfnvmrddlon = 'vmrddlon.'//trim(suffix2t09)//'.b'
853 
854 c - Due to a bug in GMT, the call will use "-I*m" rather
855 c - than "-I*s". As such, convert igridsec to xgridmin.
856  xgridmin = dble(igridsec)/60.d0
857 
858 c - Do NOT generate a call to surface if there is no data
859 c - to grid.
860 
861 c - Call to grid the RMS'd latitudes
862  if(nrmslat.ne.0)then
863  write(99,501)trim(gfnssrddlat),glomn,glomx,glamn,glamx,
864  * xgridmin,trim(rfnvsrddlat),cmidlat
865  write(99,501)trim(gfnsmrddlat),glomn,glomx,glamn,glamx,
866  * xgridmin,trim(rfnvmrddlat),cmidlat
867  endif
868 
869 c - Call to grid the RMS'd longitudes
870  if(nrmslon.ne.0)then
871  write(99,501)trim(gfnssrddlon),glomn,glomx,glamn,glamx,
872  * xgridmin,trim(rfnvsrddlon),cmidlat
873  write(99,501)trim(gfnsmrddlon),glomn,glomx,glamn,glamx,
874  * xgridmin,trim(rfnvmrddlon),cmidlat
875  endif
876 
877 c - Call to grid the RMS'd ellipsoid heights
878  if(nrmseht.ne.0)
879  *write(99,501)trim(gfnsmrddeht),glomn,glomx,glamn,glamx,
880  *xgridmin,trim(rfnvmrddeht),cmidlat
881 
882 c - Calls to convert ".grd" to ".xyz"
883  if(nrmslat.ne.0)then
884  write(99,502)trim(rfnvsrddlat),trim(zfnvsrddlat)
885  write(99,502)trim(rfnvmrddlat),trim(zfnvmrddlat)
886  endif
887 
888  if(nrmslon.ne.0)then
889  write(99,502)trim(rfnvsrddlon),trim(zfnvsrddlon)
890  write(99,502)trim(rfnvmrddlon),trim(zfnvmrddlon)
891  endif
892  if(nrmseht.ne.0)
893  *write(99,502)trim(rfnvmrddeht),trim(zfnvmrddeht)
894 
895 c - Calls to convert ".xyz" to ".b"
896  if(nrmslat.ne.0)then
897  write(99,503)trim(zfnvsrddlat),trim(bfnvsrddlat)
898  write(99,503)trim(zfnvmrddlat),trim(bfnvmrddlat)
899  endif
900  if(nrmslon.ne.0)then
901  write(99,503)trim(zfnvsrddlon),trim(bfnvsrddlon)
902  write(99,503)trim(zfnvmrddlon),trim(bfnvmrddlon)
903  endif
904  if(nrmseht.ne.0)
905  *write(99,503)trim(zfnvmrddeht),trim(bfnvmrddeht)
906 
907 c - Changed to T=0.9 on 10/28/2015
908  501 format(
909  *'surface ',a,' -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s,
910  *' -I',f0.2,'m -G',a,' -T0.9 -A',s,f6.4,' -C0.01 -V -Lld')
911 
912 c 501 format(
913 c *'surface ',a,' -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s,
914 c *' -I',f0.2,'m -G',a,' -T0.4 -A',s,f6.4,' -C0.01 -V -Lld')
915 c 501 format(
916 c *'surface ',a,' -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s,
917 c *' -I',f0.2,'m -G',a,' -T0.4 -A',s,f6.4,' -C0.01 -V')
918  502 format(
919  *'grd2xyz ',a,' -bo3f > ',a)
920  503 format(
921  *'xyz2b << !',/,
922  *a,/,a,/,'!')
923 
924 c --------------------------------------------------------------
925 c - GMT Batch file: CREATE TRANSFORMATION ERROR GRID
926 c - Added 11/9/2015
927 c
928 c - See DRU-11, p. 145-148: The total transformation error
929 c - grid (representing a "standard deviation of the transformation
930 c - at any grid location") will be a combination of two grids
931 c - that have been called (in shorthand):
932 c 1) The "method noise" grid
933 c 2) The "data noise" grid
934 c * The "method noise" grid represents the noise that comes
935 c from not knowing the perfect tension at which to pull a
936 c transformation grid sheet. It was created in batch
937 c file "gmtbat02...", as created by mymedian5.f
938 c * The "data noise" grid represents the RMS mismatch between
939 c all non-outlier vectors (both kept and dropped) and the
940 c transformation grid (as pulled at T=0.4). It was created
941 c in batch file "gmtbat05..." as created by myrms.f.
942 c The final "transformation error grid" will be a
943 c sort of convolution of the two above grids, as follows:
944 c Let "TE(i,j)" be the transformation error at some
945 c grid node "i,j", and "MN(i,j)" be the method noise
946 c at grid node i,j and "DN(i,j)" be the data noise at
947 c grid node i,j. All three should be in the same units
948 c (either meters or arcseconds). The final "TE" grid
949 c will have these values at every grid node:
950 c TE(i,j) = SQRT[ MN(i,j)**2 + DN(i,j)**2 ]
951 c --------------------------------------------------------------
952 c - The "method noise" grids
953  sadbfnvmtcdlat1000 = 'vmtcdlat.'//trim(suffix2d3)//'.b'
954  sadbfnvstcdlat1000 = 'vstcdlat.'//trim(suffix2d3)//'.b'
955  sadbfnvmtcdlon1000 = 'vmtcdlon.'//trim(suffix2d3)//'.b'
956  sadbfnvstcdlon1000 = 'vstcdlon.'//trim(suffix2d3)//'.b'
957  sadbfnvmtcdeht1000 = 'vmtcdeht.'//trim(suffix2d3)//'.b'
958 
959 c - The "total error" grids in ".b" format
960  bfnvsetelat ='vsetelat.'//trim(suffix2)//'.b'
961  bfnvsetelon ='vsetelon.'//trim(suffix2)//'.b'
962  bfnvmetelat ='vmetelat.'//trim(suffix2)//'.b'
963  bfnvmetelon ='vmetelon.'//trim(suffix2)//'.b'
964  bfnvmeteeht ='vmeteeht.'//trim(suffix2)//'.b'
965 
966 c - The "total error" grids in ".grd" format
967  rfnvsetelat ='vsetelat.'//trim(suffix2)//'.grd'
968  rfnvsetelon ='vsetelon.'//trim(suffix2)//'.grd'
969  rfnvmetelat ='vmetelat.'//trim(suffix2)//'.grd'
970  rfnvmetelon ='vmetelon.'//trim(suffix2)//'.grd'
971  rfnvmeteeht ='vmeteeht.'//trim(suffix2)//'.grd'
972 
973 
974  if(nlat.ne.0)then
975 c------------------
976 c - LATITUDE/METERS
977 c------------------
978 c - Square Method Noise:
979  write(99,801)trim(sadbfnvmtcdlat1000),
980  * trim(sadbfnvmtcdlat1000),
981  * trim(sadbfnvmtcdlat1000)
982 c - Square Data Noise:
983  write(99,802)trim(bfnvmrddlat),trim(bfnvmrddlat),
984  * trim(bfnvmrddlat)
985 c - Add the two squares together:
986  write(99,803)
987 c - Square root of the sum to give the final total transformation error grid
988  write(99,804)trim(bfnvmetelat),trim(bfnvmetelat),
989  * trim(bfnvmetelat)
990 c - Remove all dummy files
991  write(99,805)
992 c------------------
993 c - LATITUDE/ARCSECONDS
994 c------------------
995 c - Square Method Noise:
996  write(99,801)trim(sadbfnvstcdlat1000),
997  * trim(sadbfnvstcdlat1000),
998  * trim(sadbfnvstcdlat1000)
999 c - Square Data Noise:
1000  write(99,802)trim(bfnvsrddlat),trim(bfnvsrddlat),
1001  * trim(bfnvsrddlat)
1002 c - Add the two squares together:
1003  write(99,803)
1004 c - Square root of the sum to give the final total transformation error grid
1005  write(99,804)trim(bfnvsetelat),trim(bfnvsetelat),
1006  * trim(bfnvsetelat)
1007 c - Remove all dummy files
1008  write(99,805)
1009  endif
1010 
1011  if(nlon.ne.0)then
1012 c------------------
1013 c - LONGITUDE/METERS
1014 c------------------
1015 c - Square Method Noise:
1016  write(99,801)trim(sadbfnvmtcdlon1000),
1017  * trim(sadbfnvmtcdlon1000),
1018  * trim(sadbfnvmtcdlon1000)
1019 c - Square Data Noise:
1020  write(99,802)trim(bfnvmrddlon),trim(bfnvmrddlon),
1021  * trim(bfnvmrddlon)
1022 c - Add the two squares together:
1023  write(99,803)
1024 c - Square root of the sum to give the final total transformation error grid
1025  write(99,804)trim(bfnvmetelon),trim(bfnvmetelon),
1026  * trim(bfnvmetelon)
1027 c - Remove all dummy files
1028  write(99,805)
1029 c------------------
1030 c - LONGITUDE/ARCSECONDS
1031 c------------------
1032 c - Square Method Noise:
1033  write(99,801)trim(sadbfnvstcdlon1000),
1034  * trim(sadbfnvstcdlon1000),
1035  * trim(sadbfnvstcdlon1000)
1036 c - Square Data Noise:
1037  write(99,802)trim(bfnvsrddlon),trim(bfnvsrddlon),
1038  * trim(bfnvsrddlon)
1039 c - Add the two squares together:
1040  write(99,803)
1041 c - Square root of the sum to give the final total transformation error grid
1042  write(99,804)trim(bfnvsetelon),trim(bfnvsetelon),
1043  * trim(bfnvsetelon)
1044 c - Remove all dummy files
1045  write(99,805)
1046  endif
1047 
1048  if(neht.ne.0)then
1049 c------------------
1050 c - ELLIPSOID HEIGHT/METERS
1051 c------------------
1052 c - Square Method Noise:
1053  write(99,801)trim(sadbfnvmtcdeht1000),
1054  * trim(sadbfnvmtcdeht1000),
1055  * trim(sadbfnvmtcdeht1000)
1056 c - Square Data Noise:
1057  write(99,802)trim(bfnvmrddeht),trim(bfnvmrddeht),
1058  * trim(bfnvmrddeht)
1059 c - Add the two squares together:
1060  write(99,803)
1061 c - Square root of the sum to give the final total transformation error grid
1062  write(99,804)trim(bfnvmeteeht),trim(bfnvmeteeht),
1063  * trim(bfnvmeteeht)
1064 c - Remove all dummy files
1065  write(99,805)
1066  endif
1067 
1068 
1069  801 format(
1070  *'# ------------------------------',/,
1071  *'# Squaring the Method Noise Grid: ',a,' = dummy1',/,
1072  *'# ------------------------------',/,
1073  *'echo Squaring the Method Noise Grid: ',a,' = dummy1',/,
1074  *'gsqr << !',/,
1075  *a,/,'dummy1',/,'!')
1076 
1077  802 format(
1078  *'# ------------------------------',/,
1079  *'# Squaring the Data Noise Grid: ',a,' = dummy2',/,
1080  *'# ------------------------------',/,
1081  *'echo Squaring the Data Noise Grid: ',a,' = dummy2',/,
1082  *'gsqr << !',/,
1083  *a,/,'dummy2',/,'!')
1084  803 format(
1085  *'# ------------------------------',/,
1086  *'# Adding dummy1 and dummy2 = dummy3',/,
1087  *'# ------------------------------',/,
1088  *'echo Adding dummy1 and dummy2 = dummy3',/,
1089  *'addem << !',/,
1090  *'dummy1',/,'dummy2',/,'dummy3',/,'!')
1091  804 format(
1092  *'# ------------------------------',/,
1093  *'# SquareRoot dummy3 to get Trans. Error Grid: ',a,/,
1094  *'# ------------------------------',/,
1095  *'echo SquareRoot dummy3 to get Trans. Error Grid: ',a,/,
1096  *'gsqrt << !',/,
1097  *'dummy3',/,a,/,'!')
1098  805 format(
1099  *'# ------------------------------',/,
1100  *'# Removing dummy1, dummy2, dummy3',/,
1101  *'# ------------------------------',/,
1102  *'echo Removing dummy1, dummy2, dummy3',/,
1103  *'rm -f dummy1',/,
1104  *'rm -f dummy2',/,
1105  *'rm -f dummy3')
1106 
1107 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
1108 c - 2016 07 01 : Mask/MASK/mask for the HARN/FBN/CONUS situation
1109 c - See DRU-12, p. 41
1110 c - Basically: Once the "...ete...b" files are created,
1111 c - move them to be "...ete...b.premask". Then apply
1112 c - the mask to this "premask" grid (densify to 30",
1113 c - convolve with 30" mask, then decimate back to original
1114 c - grid spacing) and give it the standard "...ete...b" file name.
1115 c - In this way, when we convert the ".b" to ".grd" in the
1116 c - next section of code, it'll be for the MASKED version
1117 c - of the ".b" grid (but again, ONLY for the HARN/FBN/CONUS
1118 c - situation)
1119  if(trim(olddtm).eq.'nad83_harn' .and.
1120  * trim(newdtm).eq.'nad83_fbn' .and.
1121  * trim(region).eq.'conus' )then
1122 
1123  write(99,601)
1124 
1125  if(nlat.ne.0)then
1126  write(99,602)trim(bfnvmetelat),
1127  * trim(bfnvmetelat)//'.premask',
1128  * trim(bfnvmetelat)//'.premask',
1129  * trim(bfnvmetelat),
1130  * igridsec/30,igridsec/30
1131 
1132  write(99,602)trim(bfnvsetelat),
1133  * trim(bfnvsetelat)//'.premask',
1134  * trim(bfnvsetelat)//'.premask',
1135  * trim(bfnvsetelat),
1136  * igridsec/30,igridsec/30
1137  endif
1138 
1139  if(nlon.ne.0)then
1140  write(99,602)trim(bfnvmetelon),
1141  * trim(bfnvmetelon)//'.premask',
1142  * trim(bfnvmetelon)//'.premask',
1143  * trim(bfnvmetelon),
1144  * igridsec/30,igridsec/30
1145 
1146  write(99,602)trim(bfnvsetelon),
1147  * trim(bfnvsetelon)//'.premask',
1148  * trim(bfnvsetelon)//'.premask',
1149  * trim(bfnvsetelon),
1150  * igridsec/30,igridsec/30
1151  endif
1152 
1153  if(neht.ne.0)then
1154  write(99,602)trim(bfnvmeteeht),
1155  * trim(bfnvmeteeht)//'.premask',
1156  * trim(bfnvmeteeht)//'.premask',
1157  * trim(bfnvmeteeht),
1158  * igridsec/30,igridsec/30
1159  endif
1160 
1161  endif
1162 
1163 
1164  601 format('echo Applying MASK for HARN FBN CONUS ete grids')
1165  602 format(
1166  *'mv ',a,' ',a,/,
1167  *'regrd2 << !',/,
1168  *a,/,
1169  *'dummy1.b',/,
1170  *'3121',/,
1171  *'7081',/,
1172  *'!',/,
1173  *'convlv << !',/
1174  *'dummy1.b',/,
1175  *'Masks/mask.harnfbn.30.b',/,
1176  *'dummy2.b',/,
1177  *'!',/,
1178  *'decimate << !',/
1179  *'dummy2.b',/,
1180  *a,/,
1181  *i3,/,
1182  *i3,/,
1183  *'!',/,
1184  *'rm -f dummy1.b',/,
1185  *'rm -f dummy2.b')
1186 
1187 c --------------------------------------------------------------
1188 c - GMT Batch file: Convert transformation error grid to "grd"
1189 c - format so I can color plot it later.
1190 c --------------------------------------------------------------
1191 
1192  if(nlat.ne.0)then
1193  write(99,507)trim(bfnvmetelat)
1194  * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvmetelat)
1195  write(99,507)trim(bfnvsetelat)
1196  * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvsetelat)
1197  endif
1198 
1199  if(nlon.ne.0)then
1200  write(99,507)trim(bfnvmetelon)
1201  * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvmetelon)
1202  write(99,507)trim(bfnvsetelon)
1203  * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvsetelon)
1204  endif
1205 
1206  if(neht.ne.0)then
1207  write(99,507)trim(bfnvmeteeht)
1208  * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvmeteeht)
1209 
1210  endif
1211  507 format(
1212  *'b2xyz << !',/,a,/,'!',/,
1213  *'xyz2grd temp.xyz -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s,
1214  *' -I',f0.2,'m -bi3f -G',a,/,
1215  *'rm -f temp.xyz')
1216 
1217 
1218 
1219 
1220 
1221 
1222 
1223 
1224 
1225 
1226 
1227 
1228 
1229 
1230 
1231  write(99,1031)trim(gmtfile)
1232  1031 format('echo END batch file ',a)
1233  close(99)
1234 
1235  write(6,9999)
1236  9999 format('END program myrms.f')
1237 
1238 
1239 
1240  2006 format('Cell ',i8,' Poscode: ',i12)
1241  2007 format(6x,'Pt : ',i5)
1242  end
1243 c
1244 c -----------------------------------------------
1245 c
1246  subroutine getrms(zs,max,inumlo,inumhi,indx,
1247  *maxppcell,xlat,xlon,alat,alon,rval)
1249 c - As of 10/28/2015, no longer using average
1250 c - lat/lon, but rather center of cell, to
1251 c - register the RMS vector
1252 
1253  real*8 xlat(max),xlon(max),zs(max)
1254  integer*4 indx(max),indx2(maxppcell)
1255  integer*4 inumlo,inumhi
1256 
1257  real*8 alat,alon,rval
1258 
1259 c write(6,1)inumlo,inumhi
1260  1 format('Inside getrms',/,
1261  *'inumlo,inumhi = ',i8,1x,i8)
1262 
1263 c - When the whole set of "zs" data is sorted by
1264 c - "whichcell", then the "indx" values are our
1265 c - key to that sort. The "zs" data remains
1266 c - in "read order" (1-nkt), but the indx values
1267 c - (also in "read order" point to the sorted
1268 c - order.
1269 
1270 c - Because we are NOT attempting any kind of
1271 c - median filter, there is no need to sort
1272 c - the "mini vector" of data in this cell.
1273 c - Just compute the average lat, average lon
1274 c - and RMS value and return.
1275 c - Change as of 10/28/2015: Do not
1276 c - use the avelat/avelon. Instead,
1277 c - register the RMS vector to the center of
1278 c - the cell. (See DRU-11, p. 145)
1279 
1280  rval = 0.d0
1281  alat = 0.d0
1282  alon = 0.d0
1283 
1284 
1285  nval = inumhi-inumlo+1
1286  iq = 0
1287  do 4 i=inumlo,inumhi
1288  iq = iq + 1
1289  rval = rval + zs(indx(i))**2
1290  alat = alat + xlat(indx(i))
1291  alon = alon + xlon(indx(i))
1292  4 continue
1293  rval = dsqrt(rval/nval)
1294  alat = alat/nval
1295  alon = alon/nval
1296 
1297  return
1298  end
1299 
1300  include 'Subs/getgridbounds.f'
1301  include 'Subs/indexxi.for'
1302  include 'Subs/onzd2.f'
1303 
subroutine indexxi(n, nd, arr, indx)
Subroutine to perform ?? indexing on integer data.
Definition: indexxi.for:20
real *8 function onzd2(x)
Function to round a digit to one significant figure (one non zero digit), double precision.
Definition: onzd2.f:32
subroutine getrms(zs, max, inumlo, inumhi, indx, maxppcell, xlat, xlon, alat, alon, rval)
Definition: myrms.f:1248
subroutine getgridbounds(region, xn, xs, xw, xe)
Subroutine to collect up the GRID boundaries for use in creating NADCON 5.
Definition: getgridbounds.f:26
program myrms
Part of the NADCON5 build process, generates gmtbat05
Definition: myrms.f:97