program checkgrid c - 2019 05 20: Added "gistat" as a complement to "stat" c - See DRU-13, p. 4-6 c - 2016 05 18 : Converted from NADCON 5 to VERTCON 3 c-------------------------------------------------------- c - 2016 08 26: Changed "getmapbounds" to bring in a better way c of computing the reference vector location and added a new c variable for its label c - Also changing the call to "getmapbounds" to give it "olddatum" and "newdatum" c - to aide in filtering out things like the Saint regions in Alaska c - for unsupported transformations. c - 2016 07 29: Scrapped personal placement of vectors and just let them c - sit outside/below the map c - 2016 07 21: c - Added code to allow for optional placement of reference vectors, coming from c "map.parameters" as read in subroutine "getmapbounds" c - 2015 10 08 - Added HOR output in M and S for Lat/Lon to both "gi" and "dd" vector output. c - Combined: v(m/s)(a/t/d)(gi/dd)lat... c - v(m/s)(a/t/d)(gi/dd)lon... c - into : v(m/s)(a/t/d)(gi/dd)hor... c c c c - Program begun 9/10/2015 c - For use in creating NADCON5 c - Built by Dru Smith c - 1) Compare grids of dlat, dlon and doht c - to vectors of dlat, dlon and doht. c - 2) Spit out interpolated (from grid) vectors c - 3) Spit out differential (interpolated minus original) vectors. c - 4) Create a GMT batch file to plot said vectors. c - The input vectors: c - Represent *all* (outlier removed) c - vectors of dlat/dlon/doht for the c - olddatum/newdatum/region combination c - However, for the sake of understanding, the c - vectors will be read in from their "thinned" c - and "dropped" files, so that we can generate c - statistics of: c - thinned-versus-gridded c - dropped-versus-gridded c - all-versus-gridded c - This is important, since ONLY the thinned c - vectors went into the grid and it seems c - that their statistics should be better against c - the grid than the dropped vectors. Additionally, c - one might argue that the only independent check c - on the grid is the dropped agreement. We'll see. c - The input grids: c - dlat/dlon/doht grids based on thinning c - all of the vectors (see above) using c - a median thinning at some block spacing c - in arcseconds, and gridding to that same c - block spacing. c - Input to this program: c - olddtm c - newdtm c - region c - agridsec implicit real*8(a-h,o-z) parameter(maxlat=5000,maxlon=5000,maxpts=500000) parameter(maxplots=60) character*10 olddtm,newdtm,region character*3 filter c - Grid spacing, converted to characters for file names c - Presumes: c - 1) igridsec is an integer c - 2) 1 <= igridsec <= 99999 integer*4 igridsec,igridsec2 character*5 agridsec,agridsec2 character*1 mapflag character*200 suffix1,suffix2,suffix3 character*200 suffix2t04 character*200 gmtfile character*3 ele(7),el0(7) c - GMT stuff real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots) real*4 jm(maxplots) real*4 b1(maxplots),b2(maxplots) character*10 fn(maxplots) character*200 gfn real*8 lorvopc real*8 lorvog(7),lorvoggi(7) real*8 scaledd(7),scalegi(7) real*8 basedd(7),basegi(7) real*8 lorvogprime,lorvogmeters real*8 lorvogprimegi,lorvogmetersgi real*8 lorvoghorgis,lorvoghorgim real*8 lorvoghordds,lorvoghorddm c - 2016 07 21: logical lrv(maxplots) real*8 rv0x(maxplots),rv0y(maxplots),rl0y(maxplots) c - Input original THINNED vectors (outliers must be c - removed) character*200 gfnvmtcdlat,gfnvstcdlat character*200 gfnvmtcdlon,gfnvstcdlon character*200 gfnvmtcdoht character*200 gfnvmtcdhor,gfnvstcdhor c - Input original DROPPED vectors (outliers must be c - removed) character*200 gfnvmdcdlat,gfnvsdcdlat character*200 gfnvmdcdlon,gfnvsdcdlon character*200 gfnvmdcdoht character*200 gfnvmdcdhor,gfnvsdcdhor c - Input ".b" files created from the thinned vectors character*200 bfnvmtcdlat,bfnvstcdlat character*200 bfnvmtcdlon,bfnvstcdlon character*200 bfnvmtcdoht c - Logical arguments to tell me if there are no values in files logical*1 nothinlat,nothinlon,nothinoht c - Output grid-interpolated coordinate difference vector files: c - For the sake of completeness, we'll create grid-interpolated vectors c - for all of these combinations in meters (thus 15 files): c - thin, drop, all c - by c - lat(m),lat(s),lon(m),lon(s),oht(m) character*200 gfnvmagilat,gfnvmtgilat,gfnvmdgilat character*200 gfnvsagilat,gfnvstgilat,gfnvsdgilat character*200 gfnvmagilon,gfnvmtgilon,gfnvmdgilon character*200 gfnvsagilon,gfnvstgilon,gfnvsdgilon character*200 gfnvmagioht,gfnvmtgioht,gfnvmdgioht c - Added 2015 10 08 character*200 gfnvmagihor,gfnvmtgihor,gfnvmdgihor character*200 gfnvsagihor,gfnvstgihor,gfnvsdgihor c - Output double differenced (differences of coordinate differences, c - aka DDlat, DDlon or DDoht) GMT-ready files: c c - For the sake of completeness, we'll create differential vectors c - for all of these combinations in meters (thus 15 files): c - thin, drop, all c - by c - lat(m),lat(s),lon(m),lon(s),oht(m) c - c - Stats will be generated for thin, drop or all. c - We will create a batch file to plot thin, drop or all. c - Naming scheme is: c - 1-2: "vm" for "vectors in meters" c - 3: "a","t","d" or "r" for All, Thinned, Dropped or RMS'd c - DRU: Check that this works in "v2special" c - 4-5: "cd" or "gi" or "dd" for Coordinate Differences, Grid-Interpolated coordinate differences, or Double Differences c - 6-8: "oht" c - Thus, by way of example: c - vmacdoht = Vectors in meters, All Coordinate Differences, Orthometric Height c - vmagioht = Vectors in meters, All Grid-Interpolated coordinate differences, Orthometric Height c - and their difference would be: c - vmaddoht = Vectors in meters, Double Differences (Grid minus original), Orthometric Height c c - Other examples, in general: c - vmrddoht = Vectors in meters, RMS'd Double Differences, Orthometric Height (not used in this program) c - c - On the use of RMS: The RMS'd values only are applied to Double Differences c - (that is, Differences of Coordinate Differences). As such, we should c - not see an "r" in the 3rd position of the file name without "dd" in c - positions 4-5. These Double Differenes represent the difference c - between a Coordinate Difference interpolated off the GRID and a Coordinate c - Difference directly computed from input data (aka "vectors"). c - The RMS'd data is compiled from ALL vectors compared against the c - GRID, not just thinned vectors (which were used to MAKE the grid). c - Doing it only on "thinned" causes overly optimistic statistics, while c - doing it only on "dropped" ignores the fact that the thinned c - and dropped vectors are both in the public domain and used by c - the public.) c - Note, the "dropped" and "thinned" refer to those categories established c - way back in "mymedian5.f", under doit3.bat. c - Output double differenced vector files (all/thinned/dropped): character*200 gfnvmaddlat,gfnvmtddlat,gfnvmdddlat character*200 gfnvsaddlat,gfnvstddlat,gfnvsdddlat character*200 gfnvmaddlon,gfnvmtddlon,gfnvmdddlon character*200 gfnvsaddlon,gfnvstddlon,gfnvsdddlon character*200 gfnvmaddoht,gfnvmtddoht,gfnvmdddoht c - Added 2015 10 08 character*200 gfnvmaddhor,gfnvmtddhor,gfnvmdddhor character*200 gfnvsaddhor,gfnvstddhor,gfnvsdddhor c - Stats file name character*200 sfn c - The gridded data real*8 glamn,glomn,dla,dlo integer*4 nla,nlo,ikind real*4 datlats(maxlat,maxlon) real*4 datlons(maxlat,maxlon) real*4 datohtm(maxlat,maxlon) c - The vector data character*80 card c - All statistical data. c The "stat" and "kstat" cover "dd" data. c i,j,k,l: real*8 stat(5,5,3,3) c j,k: integer*4 kstat(5,3) c j: character*20 units(5) c - Fix from NADCON 5.01 (See DRU-13, p.5) real*8 gistat(5,5) c - stat(i,j,k,l) : c - kstat(j,k) : c c i = 1 => Average c i = 2 => RMS c i = 3 => Min c i = 4 => Max c i = 5 => Standard Deviation c c j = 1 => LAT, arcseconds c j = 2 => LON, arcseconds c j = 3 => OHT, meters c j = 4 => LAT, meters c j = 4 => LON, meters c c k = 1 => Thinned vectors c k = 2 => Dropped vectors c k = 3 => All vectors c c l = 1 => Bilinear Interpolation c l = 2 => Biquadratic Interpolation c l = 3 => Not currently used c l: real*8 vgrd(3),vgrdm(3),dd(3),ddm(3) character*4 type(3) c - Note that "dd" means arcsec/arcsec/meters for lat/lon/oht c - while "ddm" means meters/meters/meters for lat/lon/oht character*20 coord(3) real*8 fac(3,3) c - The output data real*8 xla0(maxpts),xlo0(maxpts) real*8 dd0(maxpts),ddm0(maxpts),gi0(maxpts),gim0(maxpts) character*6 pid0(maxpts) logical v2special c -------------------------------------- c -------------------------------------- c -------------------------------------- c -------BEGIN GUTS OF PROGRAM---------- c -------------------------------------- c -------------------------------------- c -------------------------------------- write(6,1001) 1001 format('BEGIN program checkgrid.f') c ------------------------------------------------------------------ c - Some required constants. c ------------------------------------------------------------------ c - Length of Reference Vector on Paper, in CM (lorvopc). Depends c - on having the "MEASURE_UNITS" in .gmtdefaults set to "cm" lorvopc = 1.d0 csm = 3.d0 pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 re = 6371000.d0 s2m = (1.d0/3600.d0)*d2r*re MultiplierLorvog = 2 coord(1) = 'LATITUDE' coord(2) = 'LONGITUDE' c coord(3) = 'ELLIPSOID HEIGHT' coord(3) = 'ORTHOMETRIC HEIGHT' units(1) = 'arcseconds' units(2) = 'arcseconds' units(3) = 'meters' units(4) = 'meters' units(5) = 'meters' c ------------------------------------------------------------------ c - Read in arguments from batch file c ------------------------------------------------------------------ read(5,'(a)')olddtm read(5,'(a)')newdtm read(5,'(a)')region read(5,'(a)')filter read(5,'(a)')agridsec read(5,'(a)')agridsec2 read(5,'(a)')mapflag v2special = .false. if(trim(region).eq.'conus' .and. * trim(olddtm).eq.'ngvd29' .and. * trim(newdtm).eq.'navd88') v2special = .true. c ------------------------------------------------------------------ c - Generate the suffixes used in all our files c ------------------------------------------------------------------ read(agridsec,*)igridsec read(agridsec2,*)igridsec2 c suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region) suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region) *//'.'//trim(filter) c suffix2=trim(suffix1)//'.'//trim(agridsec) suffix2=trim(suffix1)//'.'//trim(agridsec)//'.'//trim(agridsec2) suffix2t04=trim(suffix2)//'.04' suffix3=trim(suffix2)//'.'//trim(mapflag) c ------------------------------------------------------------------ c - Open the input, thinned, GMT-ready files c - "t" meaning "thinned". Briefly snoop the c - first line of each file. If it is empty, c - then skip using that particular component. c ------------------------------------------------------------------ c - Always maintain THIS order: c - Lat/Seconds 1 c - Lon/Seconds 2 c - Oht/Meters 3 c - (Hor/Seconds 4) c - (Hor/Meters 5) c - Lat/Meters 6 c - Lon/Meters 7 c gfnvstcdlat = 'vstcdlat.'//trim(suffix2) c open(21,file=gfnvstcdlat,status='old',form='formatted') c write(6,1012)trim(gfnvstcdlat) c gfnvstcdlon = 'vstcdlon.'//trim(suffix2) c open(22,file=gfnvstcdlon,status='old',form='formatted') c write(6,1012)trim(gfnvstcdlon) gfnvmtcdoht = 'vmtcdoht.'//trim(suffix2) open(23,file=gfnvmtcdoht,status='old',form='formatted') write(6,1012)trim(gfnvmtcdoht) c gfnvmtcdlat = 'vmtcdlat.'//trim(suffix2) c open(26,file=gfnvmtcdlat,status='old',form='formatted') c write(6,1012)trim(gfnvmtcdlat) c gfnvmtcdlon = 'vmtcdlon.'//trim(suffix2) c open(27,file=gfnvmtcdlon,status='old',form='formatted') c write(6,1012)trim(gfnvmtcdlon) c - Note: "nothinlat" means there were no thinned lat vectors, c - and thus no ".b" grid made from thinned vectors. As such, c - anything having to do with latitude grids is skipped. c - Same for lon or oht. c call nlines(21,nthinlat,nothinlat) c call nlines(22,nthinlon,nothinlon) call nlines(23,nthinoht,nothinoht) c - No need to do this for files 26/27 as the parallel 21 and 22 1012 format(6x,'checkgrid.f: ', *'Opening input thinned vector file ',a) c ------------------------------------------------------------------ c - Open the input, dropped, GMT-ready files c - "d" meaning "dropped" c ------------------------------------------------------------------ c gfnvsdcdlat = 'vsdcdlat.'//trim(suffix2) c open(31,file=gfnvsdcdlat,status='old',form='formatted') c write(6,1013)trim(gfnvsdcdlat) c gfnvsdcdlon = 'vsdcdlon.'//trim(suffix2) c open(32,file=gfnvsdcdlon,status='old',form='formatted') c write(6,1013)trim(gfnvsdcdlon) gfnvmdcdoht = 'vmdcdoht.'//trim(suffix2) open(33,file=gfnvmdcdoht,status='old',form='formatted') write(6,1013)trim(gfnvmdcdoht) c gfnvmdcdlat = 'vmdcdlat.'//trim(suffix2) c open(36,file=gfnvmdcdlat,status='old',form='formatted') c write(6,1013)trim(gfnvmdcdlat) c gfnvmdcdlon = 'vmdcdlon.'//trim(suffix2) c open(37,file=gfnvmdcdlon,status='old',form='formatted') c write(6,1013)trim(gfnvmdcdlon) 1013 format(6x,'checkgrid.f: ', *'Opening input dropped vector file ',a) c ------------------------------------------------------------------ c - Open the input, gridded-from-thinned .b files. If the c - thinned vector file is empty, skip trying to open the c - related ".b" file. c ------------------------------------------------------------------ if(.not.nothinlat)then c bfnvstcdlat = 'vstcdlat.'//trim(suffix2t04)//'.b' c open(11,file=bfnvstcdlat,status='old',form='unformatted') c write(6,1015)trim(bfnvstcdlat) c bfnvmtcdlat = 'vmtcdlat.'//trim(suffix2t04)//'.b' c open(16,file=bfnvmtcdlat,status='old',form='unformatted') c write(6,1015)trim(bfnvmtcdlat) endif if(.not.nothinlon)then c bfnvstcdlon = 'vstcdlon.'//trim(suffix2t04)//'.b' c open(12,file=bfnvstcdlon,status='old',form='unformatted') c write(6,1015)trim(bfnvstcdlon) c bfnvmtcdlon = 'vmtcdlon.'//trim(suffix2t04)//'.b' c open(17,file=bfnvmtcdlon,status='old',form='unformatted') c write(6,1015)trim(bfnvmtcdlon) endif if(.not.nothinoht)then bfnvmtcdoht = 'vmtcdoht.'//trim(suffix2t04)//'.b' open(13,file=bfnvmtcdoht,status='old',form='unformatted') write(6,1015)trim(bfnvmtcdoht) endif 1015 format(6x,'checkgrid.f: ', *'Opening input grid file ',a) c ------------------------------------------------------------------ c - Open the output, interpolated-from-grid coordinate difference c - vector files, GMT-ready for plotting. c - If the input thinned vector file is empty (aka there is c - no grid), skip c ------------------------------------------------------------------ if(.not.nothinlat)then c gfnvsagilat = 'vsagilat.'//trim(suffix2) c open(91,file=gfnvsagilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvsagilat) c gfnvmagilat = 'vmagilat.'//trim(suffix2) c open(96,file=gfnvmagilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvmagilat) c gfnvstgilat = 'vstgilat.'//trim(suffix2) c open(71,file=gfnvstgilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvstgilat) c gfnvmtgilat = 'vmtgilat.'//trim(suffix2) c open(76,file=gfnvmtgilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvmtgilat) c gfnvsdgilat = 'vsdgilat.'//trim(suffix2) c open(81,file=gfnvsdgilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvsdgilat) c gfnvmdgilat = 'vmdgilat.'//trim(suffix2) c open(86,file=gfnvmdgilat,status='new',form='formatted') c write(6,1017)'lat',trim(gfnvmdgilat) endif if(.not.nothinlon)then c gfnvsagilon = 'vsagilon.'//trim(suffix2) c open(92,file=gfnvsagilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvsagilon) c gfnvmagilon = 'vmagilon.'//trim(suffix2) c open(97,file=gfnvmagilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvmagilon) c gfnvstgilon = 'vstgilon.'//trim(suffix2) c open(72,file=gfnvstgilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvstgilon) c gfnvmtgilon = 'vmtgilon.'//trim(suffix2) c open(77,file=gfnvmtgilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvmtgilon) c gfnvsdgilon = 'vsdgilon.'//trim(suffix2) c open(82,file=gfnvsdgilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvsdgilon) c gfnvmdgilon = 'vmdgilon.'//trim(suffix2) c open(87,file=gfnvmdgilon,status='new',form='formatted') c write(6,1017)'lon',trim(gfnvmdgilon) endif if(.not.nothinoht)then gfnvmagioht = 'vmagioht.'//trim(suffix2) open(93,file=gfnvmagioht,status='new',form='formatted') write(6,1017)'oht',trim(gfnvmagioht) gfnvmtgioht = 'vmtgioht.'//trim(suffix2) open(73,file=gfnvmtgioht,status='new',form='formatted') write(6,1017)'oht',trim(gfnvmtgioht) gfnvmdgioht = 'vmdgioht.'//trim(suffix2) open(83,file=gfnvmdgioht,status='new',form='formatted') write(6,1017)'oht',trim(gfnvmdgioht) endif c - Added 2015 10 08 if(.not.nothinlat .and. .not.nothinlon)then c gfnvsagihor = 'vsagihor.'//trim(suffix2) c open(94,file=gfnvsagihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvsagihor) c gfnvmagihor = 'vmagihor.'//trim(suffix2) c open(95,file=gfnvmagihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvmagihor) c gfnvstgihor = 'vstgihor.'//trim(suffix2) c open(74,file=gfnvstgihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvstgihor) c gfnvmtgihor = 'vmtgihor.'//trim(suffix2) c open(75,file=gfnvmtgihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvmtgihor) c gfnvsdgihor = 'vsdgihor.'//trim(suffix2) c open(84,file=gfnvsdgihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvsdgihor) c gfnvmdgihor = 'vmdgihor.'//trim(suffix2) c open(85,file=gfnvmdgihor,status='new',form='formatted') c write(6,1017)'hor',trim(gfnvmdgihor) endif 1017 format(6x,'checkgrid.f: ', *'Opening output grid-interpolated ',a,' vector file :',a) c ------------------------------------------------------------------ c - Open the output, differential, GMT-ready files for c - plotting. THINNED. If the input thinned vector file c - is empty, skip c ------------------------------------------------------------------ if(.not.nothinlat)then c gfnvstddlat = 'vstddlat.'//trim(suffix2) c open(41,file=gfnvstddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvstddlat) c gfnvmtddlat = 'vmtddlat.'//trim(suffix2) c open(46,file=gfnvmtddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvmtddlat) endif if(.not.nothinlon)then c gfnvstddlon = 'vstddlon.'//trim(suffix2) c open(42,file=gfnvstddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvstddlon) c gfnvmtddlon = 'vmtddlon.'//trim(suffix2) c open(47,file=gfnvmtddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvmtddlon) endif if(.not.nothinoht)then gfnvmtddoht = 'vmtddoht.'//trim(suffix2) open(43,file=gfnvmtddoht,status='new',form='formatted') write(6,1014)'oht',trim(gfnvmtddoht) endif if(.not.nothinlon .and. .not.nothinlat)then c gfnvstddhor = 'vstddhor.'//trim(suffix2) c open(44,file=gfnvstddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvstddhor) c gfnvmtddhor = 'vmtddhor.'//trim(suffix2) c open(45,file=gfnvmtddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvmtddhor) endif c ------------------------------------------------------------------ c - Open the output, differential, GMT-ready files for c - plotting. DROPPED. c ------------------------------------------------------------------ if(.not.nothinlat)then c gfnvsdddlat = 'vsdddlat.'//trim(suffix2) c open(51,file=gfnvsdddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvsdddlat) c gfnvmdddlat = 'vmdddlat.'//trim(suffix2) c open(56,file=gfnvmdddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvmdddlat) endif if(.not.nothinlon)then c gfnvsdddlon = 'vsdddlon.'//trim(suffix2) c open(52,file=gfnvsdddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvsdddlon) c gfnvmdddlon = 'vmdddlon.'//trim(suffix2) c open(57,file=gfnvmdddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvmdddlon) endif if(.not.nothinoht)then gfnvmdddoht = 'vmdddoht.'//trim(suffix2) open(53,file=gfnvmdddoht,status='new',form='formatted') write(6,1014)'oht',trim(gfnvmdddoht) endif if(.not.nothinlat .and. .not.nothinlon)then c gfnvsdddhor = 'vsdddhor.'//trim(suffix2) c open(54,file=gfnvsdddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvsdddhor) c gfnvmdddhor = 'vmdddhor.'//trim(suffix2) c open(55,file=gfnvmdddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvmdddhor) endif c ------------------------------------------------------------------ c - Open the output, differential, GMT-ready files for c - plotting. ALL. c ------------------------------------------------------------------ if(.not.nothinlat)then c gfnvsaddlat = 'vsaddlat.'//trim(suffix2) c open(61,file=gfnvsaddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvsaddlat) c gfnvmaddlat = 'vmaddlat.'//trim(suffix2) c open(66,file=gfnvmaddlat,status='new',form='formatted') c write(6,1014)'lat',trim(gfnvmaddlat) endif if(.not.nothinlon)then c gfnvsaddlon = 'vsaddlon.'//trim(suffix2) c open(62,file=gfnvsaddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvsaddlon) c gfnvmaddlon = 'vmaddlon.'//trim(suffix2) c open(67,file=gfnvmaddlon,status='new',form='formatted') c write(6,1014)'lon',trim(gfnvmaddlon) endif if(.not.nothinoht)then gfnvmaddoht = 'vmaddoht.'//trim(suffix2) open(63,file=gfnvmaddoht,status='new',form='formatted') write(6,1014)'oht',trim(gfnvmaddoht) endif if(.not.nothinlat .and. .not.nothinlon)then c gfnvsaddhor = 'vsaddhor.'//trim(suffix2) c open(64,file=gfnvsaddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvsaddhor) c gfnvmaddhor = 'vmaddhor.'//trim(suffix2) c open(65,file=gfnvmaddhor,status='new',form='formatted') c write(6,1014)'hor',trim(gfnvmaddhor) endif 1014 format(6x,'checkgrid.f: ', *'Opening output differential ',a,' vector file :',a) c ------------------------------------------------------------------ c - Each region has officially chosen boundaries (which may c - or may not agree with the MAP boundaries). Get the c - official grid boundaries here. See DRU-11, p. 126 c ------------------------------------------------------------------ call getgridbounds(region,glamx,glamn,glomn,glomx) write(6,1004)trim(region),glamn,glamx,glomn,glomx 1004 format(6x,'checkgrid.f: Region= ',a,/, * 6x,'checkgrid.f: North = ',f12.6,/, * 6x,'checkgrid.f: South = ',f12.6,/, * 6x,'checkgrid.f: West = ',f12.6,/, * 6x,'checkgrid.f: East = ',f12.6) write(6,8002) 8002 format(50('*'),/, *6x,'checkgrid.f: BEGIN STATISTICAL REPORT ', *'(grid minus true vector diffs)') c ------------------------------------------------------- c - Read the ".b" grids into RAM c - c - See DRU-11, p. 140 for the following philosophy: c c Because the final grids in lat/lon will ONLY be c in arcseconds, all statistics should be about c arcseconds, with CONVERSION to meters on the fly, c and NOT about meter vectors compared to meter grids. c ------------------------------------------------------- c do 1 i=1,3 do 1 i=3,3 c if(i.eq.1 .and. nothinlat)goto 1 c if(i.eq.2 .and. nothinlon)goto 1 if(i.eq.3 .and. nothinoht)goto 1 ifil = 10+i read(ifil)glamn,glomn,dla,dlo,nla,nlo,ikind do 2 j=1,nla c if(i.eq.1)read(ifil)(datlats(j,k),k=1,nlo) c if(i.eq.2)read(ifil)(datlons(j,k),k=1,nlo) if(i.eq.3)read(ifil)(datohtm(j,k),k=1,nlo) 2 continue 1 continue glamx = glamn + (nla-1)*dla glomx = glomn + (nlo-1)*dlo c - HACK c write(6,*) ' After loop 1...' c ------------------------------------------------------- c - Go through the vectors. Do lat first, then lon, c - then oht. Compute DDlat, DDlon and DDoht. c - Compute region-wide statistics for: c - thinned-versus-gridded c - dropped-versus-gridded c - all-versus-gridded c - Spit out the grid interpolated lat/lon/oht values, c - as well as the double differences: c - Spit out the DDlat, DDlon and DDoht values to c - output files, for either thinned, dropped or all. c ------------------------------------------------------- c - stat(i,j,k,l) : c - kstat(j,k) : c c i = 1 => Average c i = 2 => RMS c i = 3 => Min c i = 4 => Max c i = 5 => Standard Deviation c c j = 1 => LAT, arcseconds c j = 2 => LON, arcseconds c j = 3 => OHT, meters c j = 4 => LAT, meters <- Computed by converting from arcseconds c j = 4 => LON, meters <- Computed by converting from arcseconds c c k = 1 => Thinned vectors c k = 2 => Dropped vectors c k = 3 => All vectors c c l = 1 => Bilinear Interpolation c l = 2 => Biquadratic Interpolation c l = 3 => Bicubic Spline Interpolation c - Zero stats do 3002 i=1,5 c do 3003 j=1,3 c do 3003 j=1,5 do 3003 j=3,3 do 3004 k=1,3 kstat(j,k) = 0 do 3005 l=1,3 stat(i,j,k,l) = 0.d0 if(i.eq.3)stat(i,j,k,l) = 99999.d0 if(i.eq.4)stat(i,j,k,l) =-99999.d0 3005 continue 3004 continue 3003 continue 3002 continue c - HACK c write(6,*) ' After loop 3002...' type(1) = 'THIN' type(2) = 'DROP' type(3) = 'ALL ' c ------------------------------------------------------------- c ------------------------------------------------------------- c ------------------------------------------------------------- c - MAIN LOOP c ------------------------------------------------------------- c ------------------------------------------------------------- c ------------------------------------------------------------- c In the loops below, for j,k= c 1,1 = Lat(s)/Thinned c 1,2 = Lat(s)/Dropped c 2,1 = Lon(s)/Thinned c 2,2 = Lon(s)/Dropped c 3,1 = Oht(s)/Thinned c 3,2 = Oht(s)/Dropped c - j=4 and 5, being Lat(m) and Lon(m) will be converted ON THE FLY c - from j=1 and 2 values. c 4,1 = Lat(s)/Thinned c 4,2 = Lat(s)/Dropped c 5,1 = Lon(s)/Thinned c 5,2 = Lon(s)/Dropped c do 2001 j=1,3 do 2001 j=3,3 avegiprime = 0.d0 avegimeter = 0.d0 c - VERTCON3: rmsgiprime = 0.d0 rmsgimeter = 0.d0 stdgiprime = 0.d0 stdgimeter = 0.d0 c if(j.eq.1 .and. nothinlat)goto 2001 c if(j.eq.2 .and. nothinlon)goto 2001 if(j.eq.3 .and. nothinoht)goto 2001 c - Now spin over thinned, then dropped vectors itot = 0 do 2002 k=1,2 c - HACK c write(6,*) ' Inside 2001/2002, j,k = ',j,k if(k.eq.1)ifile = 20+j ! Thinned Vectors if(k.eq.2)ifile = 30+j ! Dropped Vectors iread = 0 igood = 0 101 read(ifile,'(a)',end=103)card iread = iread + 1 c - CRITICAL, pull true values of ARCSECONDS for lat/lon c - but METERS for oht c if(j.eq.1 .or. j.eq.2)read(card, 102)xlo,xla,vvec if(j.eq.3 )read(card,1102)xlo,xla,vvec c - This shouldn't happen, but in case a vector c - is not inside the grid, skip it if(xla.gt.glamx .or. xla.lt.glamn .or. * xlo.gt.glomx .or. xlo.lt.glomn)then goto 101 endif igood = igood + 1 itot = itot + 1 c - Hold the lat, lon and PID in RAM until its time to c - spit them out. xla0(itot) = xla xlo0(itot) = xlo pid0(itot) = card(74:79) c - Do the interpolations c - "vgrd()" stores sec/sec/met for lat/lon/oht c - "vgrdm()" stores meters for everything c - Bilinear to get vgrd(1) c if(j.eq.1)call bilin(datlats,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(1)) c if(j.eq.2)call bilin(datlons,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(1)) if(j.eq.3)call bilin(datohtm,glamn,glomn,dla,dlo, * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(1)) c - Biquadratic to get vgrd(2) c if(j.eq.1)call biquad(datlats,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(2)) c if(j.eq.2)call biquad(datlons,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(2)) if(j.eq.3)call biquad(datohtm,glamn,glomn,dla,dlo, * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(2)) c - Bicubic to get vgrd(3) c if(j.eq.1)call bicubic(datlats,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(3)) c if(j.eq.2)call bicubic(datlons,glamn,glomn,dla,dlo, c * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(3)) if(j.eq.3)call bicubic(datohtm,glamn,glomn,dla,dlo, * nla,nlo,maxlat,maxlon,xla,xlo,vgrd(3)) c - Differences grd minus true (lat, lon in arcseconds; oht in meters) do 808 l = 1,3 dd(l) = vgrd(l) - vvec c - Convert seconds to meters but don't lose the seconds either c if(j.eq.1)then c vgrdm(l) = vgrd(l) * s2m c ddm(l) = dd(l) * s2m c elseif(j.eq.2)then c coslat = dcos(xla*d2r) c vgrdm(l) = vgrd(l) * s2m * coslat c ddm(l) = dd(l) * s2m * coslat c else vgrdm(l) = vgrd(l) ddm(l) = dd(l) c endif c********************************************************* c - ALL *STATISTICS* in their primary units c - (arcsec for lat/lon, meters for oht), but then c - converted to meters/meters for lat/lon after the c - fact. c - will be in arcseconds or meters appropriately. c - These stats are all for DOUBLE DIFFERENCES c c - However, I do need to tally the average absolute c - values of lat(sec), lat(m), lon(sec), lon(m) and oht(m) of c - grid-interpolated vectors so I can scale their vector c - plots correctly. We'll hack that into here... c - VERTCON3 : Need both "ave" *and* "std" of the "gi" values c - since we are doing not just vectors, but color plots c - so expanded the above hack. c********************************************************* c - k=1,2 = thin/drop; k=3 = all c - j=1,2,3 for lat-sec,lon-sec,oht-met c - j=4,5 for lat-met, lon-met c - Tally average stat(1,j,k,l) = stat(1,j,k,l) + dd(l) stat(1,j,3,l) = stat(1,j,3,l) + dd(l) c if(j.lt.3)then c stat(1,j+3,k,l) = stat(1,j+3,k,l) + ddm(l) c stat(1,j+3,3,l) = stat(1,j+3,3,l) + ddm(l) c endif c - Tally RMS: stat(2,j,k,l) = stat(2,j,k,l) + dd(l)**2 stat(2,j,3,l) = stat(2,j,3,l) + dd(l)**2 c if(j.lt.3)then c stat(2,j+3,k,l) = stat(2,j+3,k,l) + ddm(l)**2 c stat(2,j+3,3,l) = stat(2,j+3,3,l) + ddm(l)**2 c endif c - Tally Min: if(dd(l).lt.stat(3,j,k,l))stat(3,j,k,l)=dd(l) if(dd(l).lt.stat(3,j,3,l))stat(3,j,3,l)=dd(l) c if(j.lt.3)then c if(ddm(l).lt.stat(3,j+3,k,l))stat(3,j+3,k,l)=ddm(l) c if(ddm(l).lt.stat(3,j+3,3,l))stat(3,j+3,3,l)=ddm(l) c endif c - Tally Max: if(dd(l).gt.stat(4,j,k,l))stat(4,j,k,l)=dd(l) if(dd(l).gt.stat(4,j,3,l))stat(4,j,3,l)=dd(l) c if(j.lt.3)then c if(ddm(l).gt.stat(4,j+3,k,l))stat(4,j+3,k,l)=ddm(l) c if(ddm(l).gt.stat(4,j+3,3,l))stat(4,j+3,3,l)=ddm(l) c endif 808 continue c - HACK / CHOICE: For now, use biquadratically interpolated c - data as the official output. Note that c - "dd0" stores ARCSECONDS for lat or lon and METERS for oht c - but "ddm0" stores METERS for everything c dd0(igood) = dd(2) c ddm0(igood) = ddm(2) c - Double Differences: dd0(itot) = dd(2) ddm0(itot) = ddm(2) c - Grid Interpolated: gi0(itot) = vgrd(2) gim0(itot) = vgrdm(2) avegiprime = avegiprime + dabs(gi0(itot)) avegimeter = avegimeter + dabs(gim0(itot)) c - NADCON 5.01, GMT5: For the lat, lon, eht plots, we go with c - color dots now, making it important to have the AVE and STD c - of the values, NOT using ABS. But for the hor plots, the c - average absolute value is still important for reference vector c - length. As such, going to LEAVE the "avegiprime" and "avegimeter" c - values UNTOUCHED, and introduce NEW values, called c - "avegiprimeNotAbs" and "avegimeterNotAbs" and c - "rmsgiprimeNotAbs" and "rmsgimeterNotAbs" and c - "stdgiprimeNotAbs" and "stdgimeterNotAbs" avegiprimeNotAbs = avegiprimeNotAbs + gi0(itot) avegimeterNotAbs = avegimeterNotAbs + gim0(itot) rmsgiprimeNotAbs = rmsgiprimeNotAbs + gi0(itot)**2 rmsgimeterNotAbs = rmsgimeterNotAbs + gim0(itot)**2 c - VERTCON3: rmsgiprime = rmsgiprime + (dabs(gi0(itot)))**2 rmsgimeter = rmsgimeter + (dabs(gim0(itot)))**2 c - Tally counts c - Either Thinned or Dropped: kstat(j,k)=kstat(j,k) + 1 c if(j.lt.3)then c kstat(j+3,k)=kstat(j+3,k) + 1 c endif c - All: kstat(j,3)=kstat(j,3) + 1 c if(j.lt.3)then c kstat(j+3,3)=kstat(j+3,3) + 1 c endif goto 101 103 continue 2002 continue c - Double check that the counts of thinned and dropped c - add up to the count for all. kthin = kstat(j,1) kdrop = kstat(j,2) kall = kstat(j,3) c - Won't do this for j=4/5 as they should parallel j=1/2 avegiprime = avegiprime / dble(kall) avegimeter = avegimeter / dble(kall) c - VERTCON3: rmsgiprime = dsqrt(rmsgiprime / dble(kall)) rmsgimeter = dsqrt(rmsgimeter / dble(kall)) c - NADCON 5.01, GMT5: avegiprimeNotAbs = avegiprimeNotAbs / dble(kall) avegimeterNotAbs = avegimeterNotAbs / dble(kall) rmsgiprimeNotAbs = dsqrt(rmsgiprimeNotAbs / dble(kall)) rmsgimeterNotAbs = dsqrt(rmsgimeterNotAbs / dble(kall)) factor = dble(kall)/dble(kall-1) stdgiprime = dsqrt(factor*(rmsgiprime**2-avegiprime**2)) stdgimeter = dsqrt(factor*(rmsgimeter**2-avegimeter**2)) stdgiprimeNotAbs = dsqrt(factor*(rmsgiprimeNotAbs**2 * -avegiprimeNotAbs**2)) stdgimeterNotAbs = dsqrt(factor*(rmsgimeterNotAbs**2 * -avegimeterNotAbs**2)) c - Bug fix: c - PUT THIS INTO VERTCON3: gistat(1,j) = avegiprimeNotAbs gistat(2,j) = rmsgiprimeNotAbs gistat(5,j) = stdgiprimeNotAbs c write(6,*) j,kthin c write(6,*) j,kdrop c write(6,*) j,kall if(kall.ne.kthin+kdrop)then write(6,*) ' Thin+Drop .DNE. All => STOP' stop endif if(itot.ne.kall)then write(6,*) ' itot .DNE. All => STOP' stop endif c - Finalize stats for k=1,3 (thin/drop/all) do 809 k=1,3 fac(j,k) = dble(kstat(j,k))/dble(kstat(j,k)-1) c if(j.lt.3)then c fac(j+3,k) = dble(kstat(j+3,k))/dble(kstat(j+3,k)-1) c endif do 810 l=1,3 stat(1,j,k,l)=stat(1,j,k,l)/kstat(j,k) stat(2,j,k,l)=dsqrt(stat(2,j,k,l)/kstat(j,k)) stat(5,j,k,l)=dsqrt(fac(j,k)* * (stat(2,j,k,l)**2 - stat(1,j,k,l)**2)) c if(j.lt.3)then c stat(1,j+3,k,l)=stat(1,j+3,k,l)/kstat(j+3,k) c stat(2,j+3,k,l)=dsqrt(stat(2,j+3,k,l)/kstat(j+3,k)) c stat(5,j+3,k,l)=dsqrt(fac(j+3,k)* c * (stat(2,j+3,k,l)**2 - stat(1,j+3,k,l)**2)) c endif 810 continue 809 continue c ----------------------- c - At this point we've finished collecting statistics for c - one type of coordinate (lat, lon or oht), and can spit c - out the data to files. c - c - OUTPUT the "vgrd" values to the appropriate VECTOR FILE c - OUTPUT the "vgrd-vvec" values to the appropriate VECTOR FILE c - c - Use the RMS value to compute scales needed for our output vector file. c - Why not AVE? Because the AVE is expected to be VERY close to zero c - but RMS better represents overall magnitude of disagreement. c ----------------------- c c - 2,j,3,2 = rms , lat-sec/lon-sec/oht-met/lat-met/lon-met , all , biquad, double differenced c Convert the RMS to cm to go along with all over vector plots heretofore rms0prime = stat(2,j,3,2) lorvogprime = onzd2(MultiplierLorvog*rms0prime) gprime2pc = lorvopc / lorvogprime scaledd(j) = gprime2pc basedd(j) = rms0prime c if(j.lt.3)then c rms0meters = stat(2,j+3,3,2) c lorvogmeters = onzd2(MultiplierLorvog*rms0meters) c gmeters2pc = lorvopc / lorvogmeters c scaledd(j+3) = gmeters2pc c basedd(j+3) = rms0meters c endif lorvog(j) = lorvogprime c if(j.lt.3)lorvog(j+3) = lorvogmeters c - for the grid interpolated guys... lorvogprimegi = onzd2(MultiplierLorvog*avegiprime) gprime2pcgi = lorvopc / lorvogprimegi scalegi(j) = gprime2pcgi basegi(j) = avegiprime c if(j.lt.3)then c lorvogmetersgi = onzd2(MultiplierLorvog*avegimeter) c gmeters2pcgi = lorvopc / lorvogmetersgi c scalegi(j+3) = gmeters2pcgi c basegi(j+3) = avegimeter c endif lorvoggi(j) = lorvogprimegi c if(j.lt.3)lorvoggi(j+3) = lorvogmetersgi c ------------------------------------------------------------------ c - pvlon and pvlat are the percentage of the total lon/lat c - span of the plot, from which the reference vector c - begins, starting with the Lower Left corner c ------------------------------------------------------------------ pvlon = 10.d0 pvlat = 10.d0 pvlat = (pvlat / 100.d0) pvlon = (pvlon / 100.d0) c - Set the output file number ifthin = 40+j ifdrop = 50+j ifall = 60+j c - Spin over all "grid minus vector" values, and spit the differences c - into fresh vector files. Note that THIN data is first, then DROPPED do 813 ipt=1,kall if(ipt.le.kthin)then ifout1 = ifthin else ifout1 = ifdrop endif ifout2 = ifall c - Get the Azimuth for this point, double differenced version: c - NADCON 5.01: With GMT5 we now do color-dots, NOT directional vectors, c - so the need for an azimuth is gone unless we are dealing with "hor" c - data, but we are NOT c if(j.eq.1 .or. j.eq.3)then c if(dd0(ipt).lt.0)then c az=180.d0 c else c az=0.d0 c endif c else c if(dd0(ipt).lt.0)then c az=270.d0 c else c az=90.d0 c endif c endif c - Get the Azimuth for this point, grid interpolated version: c if(j.eq.1 .or. j.eq.3)then c if(gi0(ipt).lt.0)then c azgi=180.d0 c else c azgi=0.d0 c endif c else c if(gi0(ipt).lt.0)then c azgi=270.d0 c else c azgi=90.d0 c endif c endif c - Get the scaled value for this point c - Because we aren't looping over j through 4/5, we catch c - the j=4/5 values while hitting j=1/2: vcprime = dabs(dd0(ipt)*gprime2pc) c if(j.lt.3)then c vcmeters = dabs(ddm0(ipt) * gmeters2pc) c endif giprime = dabs(gi0(ipt)*gprime2pcgi) c if(j.lt.3)then c gimeters = dabs(gim0(ipt) * gmeters2pcgi) c endif c ------------------------ c - Spit out fresh "grid interpolated" as well as "double difference" vectors c c - NADCON 5.01/GMT5: We switched how our calls to "psxy" work when moving to GMT5. Now c - just spit out ONLY lon/lat/val/PID, but be careful...put arcseconds in the c - arcseconds area and meters in the meters area, since later calls might be c - reading the exact area of a vector "card" looking for that value in that spot... c - " 140" will now be the format ONLY for METERS write-out c - "1140" will now be the format ONLY for ARCSECONDS write-out c - "2140" will be the old "140", full vector format for "hor" data c ------------------------ if(j.eq.1 .or. j.eq.2)then c - Write out for j=1,2 to thinned/dropped in primary units, double differenced: write(ifout1 ,140)xlo0(ipt),xla0(ipt),az,vcprime, * dd0(ipt),ddm0(ipt),pid0(ipt) c - Write out for j=1,2 to thinned/dropped in primary units, grid-interpolated: write(ifout1+30,140)xlo0(ipt),xla0(ipt),azgi,giprime, * gi0(ipt),gim0(ipt),pid0(ipt) c - Write out for j=1,2 to all in primary units, double differenced: write(ifout2 ,140)xlo0(ipt),xla0(ipt),az,vcprime, * dd0(ipt),ddm0(ipt),pid0(ipt) c - Write out for j=1,2 to all in primary units, grid-interpolated: write(ifout2+30,140)xlo0(ipt),xla0(ipt),azgi,giprime, * gi0(ipt),gim0(ipt),pid0(ipt) c - Write out for j=4,5 to thinned/dropped in meters, double differenced: write(ifout1+5 ,140)xlo0(ipt),xla0(ipt),az,vcmeters, * dd0(ipt),ddm0(ipt),pid0(ipt) c - Write out for j=4,5 to thinned/dropped in meters, grid interpolated: write(ifout1+35,140)xlo0(ipt),xla0(ipt),azgi,gimeters, * gi0(ipt),gim0(ipt),pid0(ipt) c - Write out for j=4,5 to all in meters, double differenced: write(ifout2+5 ,140)xlo0(ipt),xla0(ipt),az,vcmeters, * dd0(ipt),ddm0(ipt),pid0(ipt) c - Write out for j=4,5 to all in meters, grid interpolated: write(ifout2+35,140)xlo0(ipt),xla0(ipt),azgi,gimeters, * gi0(ipt),gim0(ipt),pid0(ipt) endif c - VERTCON3: We switched how our calls to "psxy" work when moving to GMT5. Now c - just spit out ONLY lon/lat/val/PID if(j.eq.3)then c - Write out for j=3 to thinned/dropped in primary units, double differenced: c write(ifout1 ,140)xlo0(ipt),xla0(ipt),az,vcprime, c * 0.d0 ,ddm0(ipt),pid0(ipt) write(ifout1 ,140)xlo0(ipt),xla0(ipt), * ddm0(ipt),pid0(ipt) c - Write out for j=3 to thinned/dropped in primary units, grid interpolated: c write(ifout1+30,140)xlo0(ipt),xla0(ipt),azgi,giprime, c * 0.d0 ,gim0(ipt),pid0(ipt) write(ifout1+30,140)xlo0(ipt),xla0(ipt), * gim0(ipt),pid0(ipt) c - Write out for j=3 to all in primary units, double differenced: c write(ifout2 ,140)xlo0(ipt),xla0(ipt),az,vcprime, c * 0.d0 ,ddm0(ipt),pid0(ipt) write(ifout2 ,140)xlo0(ipt),xla0(ipt), * ddm0(ipt),pid0(ipt) c - Write out for j=3 to all in primary units, grid interpolated: c write(ifout2+30,140)xlo0(ipt),xla0(ipt),azgi,giprime, c * 0.d0 ,gim0(ipt),pid0(ipt) write(ifout2+30,140)xlo0(ipt),xla0(ipt), * gim0(ipt),pid0(ipt) endif 813 continue c 140 format(f16.10,1x,f15.10,1x,f6.2,1x,f12.2,1x,f5.1,1x,f9.5,1x,a6) c - VERTCON3: We switched how our calls to "psxy" work when moving to GMT5. Now c - just spit out ONLY lon/lat/val/PID c 140 format(f16.10,1x,f15.10,1x,f6.2,1x,f12.2,1x,f9.5,1x,f9.3,1x,a6) 140 format(f16.10,1x,f15.10,1x,6x ,1x,12x ,1x,9x ,1x,f9.3,1x,a6) c ------------------------------------------------------------------ c - Write out a clear report. c ------------------------------------------------------------------ write(6,8000)trim(coord(j)),units(j) do 2003 k=1,3 write(6,8001)type(k), * (kstat(j,k),l=1,3), * (stat(1,j,k,l),l=1,3), * (stat(2,j,k,l),l=1,3), * (stat(5,j,k,l),l=1,3), * (stat(3,j,k,l),l=1,3), * (stat(4,j,k,l),l=1,3) 2003 continue c if(j.lt.3)then c write(6,8000)trim(coord(j)),units(j+3) c do 2004 k=1,3 c write(6,8001)type(k), c * (kstat(j+3,k),l=1,3), c * (stat(1,j+3,k,l),l=1,3), c * (stat(2,j+3,k,l),l=1,3), c * (stat(5,j+3,k,l),l=1,3), c * (stat(3,j+3,k,l),l=1,3), c * (stat(4,j+3,k,l),l=1,3) c2004 continue c endif 2001 continue c - HACK c write(6,*) ' After loop 2001...' c 102 format(f16.10,1x,f15.10,1x, 6x,1x, 12x,1x, 5x,1x,f9.5,1x,6x) c - To pull arcseconds: 102 format(f16.10,1x,f15.10,1x, 6x,1x, 12x,1x,f9.5,1x, 9x,1x,a6) c - To pull meters: 1102 format(f16.10,1x,f15.10,1x, 6x,1x, 12x,1x, 9x,1x,f9.3,1x,a6) 8000 format( *//,25x,a,' in ',a,/, *'Type',1x,'Stat', *7x,'Bilinear',5x,'Biquadratic', *9x,'Bicubic ') 8001 format(a4,1x, *'Num:',3(i15 ,1x),/, *5x,'Ave:',3(f15.6,1x),/, *5x,'RMS:',3(f15.6,1x),/, *5x,'STD:',3(f15.6,1x),/, *5x,'MIN:',3(f15.6,1x),/, *5x,'MAX:',3(f15.6,1x),/, *50('-')) write(6,8003) 8003 format( *6x,'checkgrid.f: END STATISTICAL REPORT ', *'(grid minus true vector diffs)',/, *50('*')) c ----------------------------------------------- c - CREATE THE HORIZONTAL VECTOR FILES c - 2015 10 08 c ----------------------------------------------- write(6,8004) 8004 format( *6x,'checkgrid.f: Generating HORIZONTAL vector files', *' from lat/lon vector files') 5072 format('FATAL in checkgrid: Mismatched PID for hor') 5073 format('FATAL in checkgrid: Mismatched LON for hor') 5074 format('FATAL in checkgrid: Mismatched LAT for hor') do 5080 i=1,5 if(i.eq.3)goto 5080 c write(6,*)i,'basedd() = ',basedd(i) 5080 continue do 5081 i=1,5 if(i.eq.3)goto 5081 c write(6,*)i,'basegi() = ',basegi(i) 5081 continue baseddhors = sqrt(basedd(1)**2+basedd(2)**2) baseddhorm = sqrt(basedd(4)**2+basedd(5)**2) basegihors = sqrt(basegi(1)**2+basegi(2)**2) basegihorm = sqrt(basegi(4)**2+basegi(5)**2) lorvoghorddm = onzd2(MultiplierLorvog*baseddhorm) gm2pchordd = lorvopc / lorvoghorddm lorvoghordds = onzd2(MultiplierLorvog*baseddhors) gs2pchordd = lorvopc / lorvoghordds lorvoghorgim = onzd2(MultiplierLorvog*basegihorm) gm2pchorgi = lorvopc / lorvoghorgim lorvoghorgis = onzd2(MultiplierLorvog*basegihors) gs2pchorgi = lorvopc / lorvoghorgis lorvog(6) = lorvoghordds lorvog(7) = lorvoghorddm lorvoggi(6) = lorvoghorgis lorvoggi(7) = lorvoghorgim convgs2pcddlat = lorvopc / lorvog(1) convgs2pcddlon = lorvopc / lorvog(2) convgm2pcddlat = lorvopc / lorvog(4) convgm2pcddlon = lorvopc / lorvog(5) convgs2pcgilat = lorvopc / lorvoggi(1) convgs2pcgilon = lorvopc / lorvoggi(2) convgm2pcgilat = lorvopc / lorvoggi(4) convgm2pcgilon = lorvopc / lorvoggi(5) convgs2pcddlat = scaledd(1) convgs2pcddlon = scaledd(2) convgm2pcddlat = scaledd(4) convgm2pcddlon = scaledd(5) convgs2pcgilat = scalegi(1) convgs2pcgilon = scalegi(2) convgm2pcgilat = scalegi(4) convgm2pcgilon = scalegi(5) convgs2pcddhor = dsqrt(convgs2pcddlat**2 + convgs2pcddlon**2) convgm2pcddhor = dsqrt(convgm2pcddlat**2 + convgm2pcddlon**2) convgs2pcgihor = dsqrt(convgs2pcgilat**2 + convgs2pcgilon**2) convgm2pcgihor = dsqrt(convgm2pcgilat**2 + convgm2pcgilon**2) convgs2pcddhor = gs2pchordd convgm2pcddhor = gm2pchordd convgs2pcgihor = gs2pchorgi convgm2pcgihor = gm2pchorgi c - VERTCON3: goto 8071 do 5070 i=40,90,10 do 5071 j=1,2 if(i.le.60 .and. j.eq.1)conv = convgs2pcddhor if(i.le.60 .and. j.eq.2)conv = convgm2pcddhor if(i.ge.70 .and. j.eq.1)conv = convgs2pcgihor if(i.ge.70 .and. j.eq.2)conv = convgm2pcgihor ilat = i+1 + (j-1)*5 ilon = i+2 + (j-1)*5 ihor = i+4 + (j-1)*1 c write(6,*) ilat,ilon,ihor rewind(ilat) rewind(ilon) 5078 read(ilat,140,end=5071)ylo1,yla1,az1,vc1,sec1,xmet1,pid1 read(ilon,140)ylo2,yla2,az2,vc2,sec2,xmet2,pid2 if(pid1.ne.pid2)then write(6,5072) stop endif if(ylo1.ne.ylo2)then write(6,5073) stop endif if(yla1.ne.yla2)then write(6,5074) stop endif azhor = datan2(xmet2,xmet1)/d2r if(azhor.lt.0.d0)azhor=azhor+360.d0 dhors = dsqrt(sec1**2 + sec2**2) dhorm = dsqrt(xmet1**2 + xmet2**2) if(j.eq.1)vchor = conv * dhors if(j.eq.2)vchor = conv * dhorm write(ihor,140)ylo1,yla1,azhor,vchor,dhors,dhorm,pid1 goto 5078 5071 continue 5070 continue 8071 continue c ------------------------------------------------------------------ c - Write out statistics that I'll pick up later when c - I run "makeplotfiles03.f" -- just the "rmslatm", "rmslonm" c - and "rmsohtm" values, used to scale vector plots. c ------------------------------------------------------------------ sfn = 'dvstats.'//trim(suffix2) open(2,file=sfn,status='new',form='formatted') c 2,j,3,2 = RMS/j/All/Biquad, c - where j=1/2/3 => lat-sec/lon-sec/oht-met c - j=4/5 => lat-met/lon-met c - VERTCON3: Write ONLY the OHT value c write(2,777)kstat(1,3),stat(2,1,3,2) c write(2,777)kstat(2,3),stat(2,2,3,2) write(2,777)kstat(3,3),stat(2,3,3,2) c write(2,777)kstat(1,3),stat(2,4,3,2) c write(2,777)kstat(2,3),stat(2,5,3,2) c - 2015 10 09: rhors = sqrt(stat(2,1,3,2)**2 + stat(2,2,3,2)**2) rhorm = sqrt(stat(2,4,3,2)**2 + stat(2,5,3,2)**2) c - VERTCON3, skip: c write(2,777)kstat(1,3),rhors c write(2,777)kstat(2,3),rhorm close(2) 777 format(i10,1x,f20.10) c ---------------------------------------------------------- c ---------------------------------------------------------- c ---------------------------------------------------------- c - Create GMT file "gmtbat04..." and fill it full of c - calls to make vector plots of thinned/dropped/all c - "grid interpolated" vectors as well as c - "double differenced" vectors. c ---------------------------------------------------------- c ---------------------------------------------------------- c ---------------------------------------------------------- c ------------------------------------------------------------------ c - GMT Batch file: Open the file. The output batch file full c - of all the GMT commands needed to create c - the plots of interest c ------------------------------------------------------------------ gmtfile = 'gmtbat04.'//trim(suffix3) open(99,file=gmtfile,status='new',form='formatted') write(6,1011)trim(gmtfile) 1011 format(6x,'checkgrid.f: Creating GMT batch file ',a) write(99,1030)trim(gmtfile) 1030 format('echo BEGIN batch file ',a) c -------------------------------------------------------------- c - GMT Batch file: Determine Number of areas to plot (=nplots) c - and Boundaries of plots and other stuff. c - Report the information out. c -------------------------------------------------------------- c - 2016 08 29 call getmapbounds(mapflag,maxplots,region,nplots, *olddtm,newdtm, *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) write(6,1006)trim(region) 1006 format(6x,'checkgrid.f: Calling getmapbounds for region ',a) c -------------------------------------------------------------- c - GMT Batch file: Like "makeplotfiles02.f", I will c - spin over "nplots" as my main loop, inside of that though, c - I will spin over j,k (as above, but letting k go from 1 to 3) c - Remember: j=1,2,3 => lat,lon,oht c - k=1,2,3 => Thinned, Dropped, All c -------------------------------------------------------------- ele(1) = 'lat' ele(2) = 'lon' ele(3) = 'oht' ele(4) = 'lat' ele(5) = 'lon' ele(6) = 'hor' ele(7) = 'hor' el0(1) = 'LAT' el0(2) = 'LON' el0(3) = 'OHT' el0(4) = 'LAT' el0(5) = 'LON' el0(6) = 'HOR' el0(7) = 'HOR' 991 format( *'# ------------------------------',/, *'# Double Difference Vector Plots',/, *'# ------------------------------',/, *'echo Double Difference Vector Plots') 992 format( *'# ------------------------------',/, *'# Grid Interpolated Vector Plots',/, *'# ------------------------------',/, *'echo Grid Interpolated Vector Plots') 990 format( *'# ------------------------------',/, *'# Plots for region: ',a,', sub-region: ',a,/, *'# ------------------------------',/, *'echo Creating plots for region: ',a,', sub-region: ',a) c---------------------------------- c - Grid-Interpolated Vector Plots c---------------------------------- write(99,992) do 4011 ij=1,nplots write(99,990)trim(region),trim(fn(ij)), * trim(region),trim(fn(ij)) c - 2016 08 26: xvlon = rv0x(ij) xvlat = rv0y(ij) xllat = rl0y(ij) xllon = xvlon c do 4012 j = 1,5 c do 4012 j = 1,7 c - VERTCON3: do 4012 j = 3,3 do 4013 k = 1,3 if(j.eq.1.and.k.eq.1)gfn = gfnvstgilat if(j.eq.1.and.k.eq.2)gfn = gfnvsdgilat if(j.eq.1.and.k.eq.3)gfn = gfnvsagilat if(j.eq.2.and.k.eq.1)gfn = gfnvstgilon if(j.eq.2.and.k.eq.2)gfn = gfnvsdgilon if(j.eq.2.and.k.eq.3)gfn = gfnvsagilon if(j.eq.3.and.k.eq.1)gfn = gfnvmtgioht if(j.eq.3.and.k.eq.2)gfn = gfnvmdgioht if(j.eq.3.and.k.eq.3)gfn = gfnvmagioht if(j.eq.4.and.k.eq.1)gfn = gfnvmtgilat if(j.eq.4.and.k.eq.2)gfn = gfnvmdgilat if(j.eq.4.and.k.eq.3)gfn = gfnvmagilat if(j.eq.5.and.k.eq.1)gfn = gfnvmtgilon if(j.eq.5.and.k.eq.2)gfn = gfnvmdgilon if(j.eq.5.and.k.eq.3)gfn = gfnvmagilon if(j.eq.6.and.k.eq.1)gfn = gfnvstgihor if(j.eq.6.and.k.eq.2)gfn = gfnvsdgihor if(j.eq.6.and.k.eq.3)gfn = gfnvsagihor if(j.eq.7.and.k.eq.1)gfn = gfnvmtgihor if(j.eq.7.and.k.eq.2)gfn = gfnvmdgihor if(j.eq.7.and.k.eq.3)gfn = gfnvmagihor if(j.le.5)then if(kstat(j,k).ne.0)then c - VERTCON3: Changed from B/W vector plots to colored dot plots c * call bwplotvc(ele(j),gfn,bw,be,bs,bn,jm,b1,b2, c * maxplots,olddtm,newdtm,region,el0(j),ij,xvlon, c * xvlat,xllon,xllat,lorvoggi(j),lorvopc,igridsec,fn) if(j.eq.3.and.k.eq.1)write(6,9010)'thin' if(j.eq.3.and.k.eq.2)write(6,9010)'dropped' if(j.eq.3.and.k.eq.3)write(6,9010)'all' 9010 format(6x,'checkgrid.f: Calling cpt for gi data: ',a) c - Bug Fix...see DRU-12, p. 5 aveZZ = gistat(1,j) stdZZ = gistat(5,j) call cpt(aveZZ,stdZZ,csm,cptloZZ, * cpthiZZ,cptinZZ) c call cpt(avegimeter,stdgimeter,csm,cptloohtmgi, c * cpthiohtmgi,cptinohtmgi) call coplotvc(ele(j),gfn,bw,be,bs,bn, * jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,el0(j),ij,xvlon, * xvlat,xllon,xllat, * lorvoggi(j),lorvopc,igridsec,igridsec2,fn, * cptloZZ,cpthiZZ,cptinZZ) c call coplotvc(ele(j),gfn,bw,be,bs,bn, c * jm,b1,b2,maxplots, c * olddtm,newdtm,region,filter,el0(j),ij,xvlon, c * xvlat,xllon,xllat, c * lorvoggi(j),lorvopc,igridsec,igridsec2,fn, c * cptloohtmgi,cpthiohtmgi,cptinohtmgi) endif else c if(kstat(1,k).ne.0 .and. kstat(2,k).ne.0) c * call bwplotvc(ele(j),gfn,bw,be,bs,bn,jm,b1,b2, c * maxplots,olddtm,newdtm,region,el0(j),ij,xvlon, c * xvlat,xllon,xllat,lorvoggi(j),lorvopc,igridsec,fn) endif 4013 continue 4012 continue 4011 continue c---------------------------------- c - Double-Difference Vector Plots c---------------------------------- write(99,991) do 3011 ij=1,nplots write(99,990)trim(region),trim(fn(ij)), * trim(region),trim(fn(ij)) c - Where is the start of the reference vector: c - 2016 07 29: Decided to scrap personalized locations c - and just put all reference vectors outside/below c - the plots c - 2016 08 26 xvlon = rv0x(ij) xvlat = rv0y(ij) xllon = xvlon xllat = rl0y(ij) c do 3012 j = 1,3 c do 3012 j = 1,5 c do 3012 j = 1,7 c - VERTCON3: do 3012 j = 3,3 do 3013 k = 1,3 if(j.eq.1.and.k.eq.1)gfn = gfnvstddlat if(j.eq.1.and.k.eq.2)gfn = gfnvsdddlat if(j.eq.1.and.k.eq.3)gfn = gfnvsaddlat if(j.eq.2.and.k.eq.1)gfn = gfnvstddlon if(j.eq.2.and.k.eq.2)gfn = gfnvsdddlon if(j.eq.2.and.k.eq.3)gfn = gfnvsaddlon if(j.eq.3.and.k.eq.1)gfn = gfnvmtddoht if(j.eq.3.and.k.eq.2)gfn = gfnvmdddoht if(j.eq.3.and.k.eq.3)gfn = gfnvmaddoht if(j.eq.4.and.k.eq.1)gfn = gfnvmtddlat if(j.eq.4.and.k.eq.2)gfn = gfnvmdddlat if(j.eq.4.and.k.eq.3)gfn = gfnvmaddlat if(j.eq.5.and.k.eq.1)gfn = gfnvmtddlon if(j.eq.5.and.k.eq.2)gfn = gfnvmdddlon if(j.eq.5.and.k.eq.3)gfn = gfnvmaddlon if(j.eq.6.and.k.eq.1)gfn = gfnvstddhor if(j.eq.6.and.k.eq.2)gfn = gfnvsdddhor if(j.eq.6.and.k.eq.3)gfn = gfnvsaddhor if(j.eq.7.and.k.eq.1)gfn = gfnvmtddhor if(j.eq.7.and.k.eq.2)gfn = gfnvmdddhor if(j.eq.7.and.k.eq.3)gfn = gfnvmaddhor if(j.le.5)then if(kstat(j,k).ne.0)then c - VERTCON3: c * call bwplotvc(ele(j),gfn,bw,be,bs,bn,jm,b1,b2, c * maxplots,olddtm,newdtm,region,el0(j),ij,xvlon, c * xvlat,xllon,xllat,lorvog(j),lorvopc,igridsec,fn) if(j.eq.3.and.k.eq.1)write(6,9011)'thin' if(j.eq.3.and.k.eq.2)write(6,9011)'dropped' if(j.eq.3.and.k.eq.3)write(6,9011)'all' 9011 format(6x,'checkgrid.f: Calling cpt for dd data: ',a) call cpt(stat(1,3,3,2),stat(5,3,3,2),csm,cptloohtmdd, * cpthiohtmdd,cptinohtmdd) call coplotvc(ele(j),gfn,bw,be,bs,bn, * jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,el0(j),ij,xvlon, * xvlat,xllon,xllat, * lorvog(j),lorvopc,igridsec,igridsec2,fn, * cptloohtmdd,cpthiohtmdd,cptinohtmdd) endif else c if(kstat(1,k).ne.0 .and. kstat(2,k).ne.0) c * call bwplotvc(ele(j),gfn,bw,be,bs,bn,jm,b1,b2, c * maxplots,olddtm,newdtm,region,el0(j),ij,xvlon, c * xvlat,xllon,xllat,lorvog(j),lorvopc,igridsec,fn) endif 3013 continue 3012 continue 3011 continue write(99,1031)trim(gmtfile) 1031 format('echo END batch file ',a) close(99) write(6,9999) 9999 format('END program checkgrid.f') end c c -------------------------------------------------------------- c subroutine nlines(ifile,num,ilogic) integer*4 ifile,num logical*1 ilogic character*1 dummy ilogic = .false. num = 0 1 read(ifile,'(a)',end=2)dummy num = num + 1 goto 1 2 if(num.eq.0)ilogic = .true. rewind(ifile) return end c c -------------------------------------------------------------- c include 'Subs/getmapbounds.f' include 'Subs/getgridbounds.f' include 'Subs/coplotvc.f' include 'Subs/cpt.f' include 'Subs/plotcoast.f' c include '/home/dru/Subs/bilin.f' include '/home/dru/Subs/biquad.f' include '/home/dru/Subs/bicubic.f' c include 'Subs/bwplotvc.f' c - VERTCON3: Don't need "onzd2" cuz it's in CPT already c include '/home/dru/Subs/onzd2.f'