program mymedian5 c - 2018 06 19: Updated to STILL thin/drop the v2special runs on "rz" (rather than "cd") data, c - BUT to also carry along a "vmtcd" and "vmdcd" file which are thinned/dropped c - identically to the "rz" data. This reflects NO change to functionality for c - non-v2special runs. c - The "all" / "cd" input unit will be 92 (input all rz = 23) c - The "thinned" / "cd" output unit will be 93 (output thin rz = 43) c - The "dropped" / "cd" output unit will be 94 (output drop rz = 63) c - 2018 05 16: Updated to allow for TWO thin/grid spacings c - to come in. In all cases except "v2special" (ngvd29/navd88/conus), c - those two spacings WILL be identical (doit3.bat will crash c - otherwise, and you won't even get here). But for "v2special" c - the two spacings are: c 1) Thin/drop/grid for RZ data c 2) The grid spacing desired of the final CD data c This allows us to grid broad regional RZ data (say 900 seconds), and c then restore to a higher res V2 grid (either the original 180 seconds or c even something else) to end up with a higher res CD grid. c To sum up: c Grid spacing 1 = RZ data c Grid spacing 2 = V2 data and CD data of final grids c c - 2018 03 26 c - Cannibalized version from NADCON 5 for creation in VERTCON 3 c c - Special code inside which will thin/grid based on "rz" rather c - than "cd" if we are in the ngvd29/navd88/conus realm... (= "v2special") c ------------------------------------------------------------------ c - 2016 09 14 c - Bug found when running Hawaii with points outside grid -- some mixup c - between "ikt" and "ipid". Changed this program as follows: It now c - REQUIRES that the incoming data has NO points outside the grid c - boundary. This has been forced by giving such points a "444" reject c - code in "makework" when creating the work file. By going through c - "makeplotfiles01" with a "444" reject, those points won't even make c - it into the "all" file, which is our input for median filtering here... c - 2016 06 29 c - Updated to insert masking commands for the "...04.b" (transformation c - grid) into gmtbat02 when working ONLY in the HARN/FBN/CONUS combination. c - 2015 10 27 c - For use in creating NADCON5 c - Built by Dru Smith c - Built from scratch, scrapping all previous c - "mymedian" programs used in making GEOCON c c c - Unlike "mymedian.f", this program is c - set up to filter/process all data c - at once in one run. Also, significant c - philosophical changes occurred, including: c - 1) Median filter on absolute horizontal length, and c - then, when we have our "kept" points, we use c - the lat and lon of those kept points to grid c - lat and lon separately. (mymedian.f actually c - sorted lat medians and lon medians separately, c - raising the very real possibility that separate c - points would go into each grid.) c - 2) Nothing RANDOM! No coin flipping, etc. It was c - viewed, for NADCON 5.0, as scientifically c - improper for the final grids to be reliant upon c - a filtering mechanism that could be different c - each time it was run. c - Program to c - 1) Run a customized block-median thinning c - algorithm on coordinate differences (horizontal c - and ellipsoid height) c - 2) Save the thinned data to GMT-ready plottable files c - 3) Save the removed data to GMT-ready plottable files c - 4) Create a GMT batch file (gmtbat02) to grid the thinned data at T=0.4 (which c - becomes the final transformation grid) and also at T=1.0 and c - T=0.0, (whose difference becomes the bases for the c - "method noise" grid) c - 5) If, and only if, we are doing the combination of c - HARN/FBN/CONUS, insert commands into gmtbat02 c - to apply a mask to the "04.b" grids (See DRU-12, p. 36-37) c - See DRU-11, p. 127 implicit real*8(a-h,o-z) c - Updated to allow up to 280,000 points to come in c - which also translates to letting 280,000 grid cells c - to be non-empty. c - VERTCON3: 430000 parameter(max=430000) c - And using a little hit or miss logic, the number c - of possible points in one cell can't be more than c - about 100, right? parameter(maxppcell=10000) c - Holding place for the values themselves, in read-order dimension zs(max) c - Pass/nopass flag for all data, in read-order logical lpass(max) c - Holding place for each record, in read-order character*80 cards(max) ! 2016 09 14 c - Index of what latititude row and longitude column c - is associated with each non-empty cell, in c - cell-found order. c - 1 <= ilapcell() <= nla c - 1 <= ilopcell() <= nlo integer*4 ilapcell(max),ilopcell(max) integer*4 poscode(max) integer*4 indx(max) c - Distance of a point from its cell's center real*8 dist(max) c - Index telling me which cell a given point is c - associated with, in read-order c - 1 <= whichcell() <= ncells integer*4 whichcell(max) c - ppcell(i) contains the count of how many values fall c - in cell "i", where "i" ratchets up from 1 to max c - every time I find a new non-empty cell. integer*4 ppcell(max) character*10 olddtm,newdtm,region character*3 filter character*200 suffix1,suffix2 character*200 suffix2t00,suffix2t04,suffix2t10 character*200 suffix2d1,suffix2d2,suffix2d3 c - Input GMT-ready files: character*200 gfncvacdoht character*200 gfnvmacdoht character*200 gfnvmarzoht character*200 gfnvmav2oht c - Output thinned GMT-ready files: character*200 gfncvtcdoht character*200 gfnvmtcdoht character*200 gfnvmtrzoht character*200 gfnvmtv2oht c - Output dropped GMT-ready files: character*200 gfncvdcdoht character*200 gfnvmdcdoht character*200 gfnvmdrzoht character*200 gfnvmdv2oht c - Output thinned GMT(surface)-ready files: character*200 gfnsmtcdoht character*200 gfnsmtrzoht character*200 gfnsmtv2oht c - To be created ".grd" files by surface character*200 rfnvmtcdoht04 character*200 rfnvmtcdoht00 character*200 rfnvmtcdoht10 character*200 rfnvmtrzoht04 character*200 rfnvmtrzoht00 character*200 rfnvmtrzoht10 c - To be created ".xyz" files by GMT's grd2xyz program character*200 zfnvmtcdoht04 character*200 zfnvmtcdoht00 character*200 zfnvmtcdoht10 character*200 zfnvmtrzoht04 character*200 zfnvmtrzoht00 character*200 zfnvmtrzoht10 c - To be created ".b" files by grd2b character*200 bfnvmtcdoht04 character*200 bfnvmtcdoht00 character*200 bfnvmtcdoht10 character*200 bfnvmtrzoht04 character*200 bfnvmtrzoht00 character*200 bfnvmtrzoht10 c - 2018 05 16: Regridded RZ file to final "igridsec2" grid for the "restore" step character*200 bfnvmtrzoht04x c - Difference grids to be used in creation of the "method noise" grid character*200 dbfnvmtcdoht1000 character*200 adbfnvmtcdoht1000 character*200 sadbfnvmtcdoht1000 character*200 sadrfnvmtcdoht1000 character*200 dbfnvmtrzoht1000 character*200 adbfnvmtrzoht1000 character*200 sadbfnvmtrzoht1000 character*200 sadrfnvmtrzoht1000 c - Dummy character*200 fused c - Output GMT-batch file: character*200 gmtfile 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 character*5 agridsec c - 2018 05 16, 2 grid spacings: integer*4 igridsec2 character*5 agridsec2 character*4 dumgridsec c - 2018 05 16: Using 2 grid spacings: character*4 dumgridsec2 character*6 pid character*80 card character*80 cardcvlat,cardcvlon,cardcvoht character*80 cardvmlat,cardvmlon,cardvmoht,cardvmhor character*80 cardvslat,cardvslon, cardvshor c - For "carrying along" the CD data while sorting/thinning/dropping on RZ data, during v2special character*80 cardvmoht0 logical v2special character*200 v2fnameb,dumfnb character*200 v2fnameg,dumfng c -------------------------------------- c -------------------------------------- c -------------------------------------- c -------BEGIN GUTS OF PROGRAM---------- c -------------------------------------- c -------------------------------------- c -------------------------------------- write(6,1001) 1001 format('BEGIN program mymedian5.f') c ------------------------------------------------------------------ c - Some required constants. c ------------------------------------------------------------------ pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 rng2std = 1.d0/1.6929d0 c ------------------------------------------------------------------ c - Read in arguments from batch file c ------------------------------------------------------------------ read(5,'(a)')olddtm read(5,'(a)')newdtm read(5,'(a)')region read(5,'(a3)')filter read(5,'(a)')agridsec read(5,'(a)')agridsec2 c ------------------------------------------------------------------ c - For NGVD29/NAVD88/CONUS only, we have a special remove/compute/restore c - process going on, with extra files, etc. Be ready to deal with that. c ------------------------------------------------------------------ v2special = .false. if(trim(olddtm).eq.'ngvd29' .and. * trim(newdtm).eq.'navd88' .and. * trim(region).eq.'conus' )v2special = .true. c ------------------------------------------------------------------ c - Generate the suffixes used in all our files c ------------------------------------------------------------------ read(agridsec,*)igridsec read(agridsec2,*)igridsec2 if(.not.v2special)then if(igridsec .ne. igridsec2)then write(6,743)igridsec,igridsec2 stop endif endif 743 format(6x,'mymedian5.f: FATAL ERROR: ', * 'Mismatched grid spacings for non-v2special run: ', * i0,1x,i0) suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region) *//'.'//trim(filter) c suffix2=trim(suffix1)//'.'//trim(agridsec) suffix2=trim(suffix1)//'.'//trim(agridsec)//'.'//trim(agridsec2) suffix2t00=trim(suffix2)//'.00' suffix2t04=trim(suffix2)//'.04' suffix2t10=trim(suffix2)//'.10' suffix2d1=trim(suffix2)//'.d1' suffix2d2=trim(suffix2)//'.d2' suffix2d3=trim(suffix2)//'.d3' c ------------------------------------------------------------------ c - Open the GMT batch file for gridding the thinned data c ------------------------------------------------------------------ gmtfile = 'gmtbat02.'//trim(suffix2) open(99,file=gmtfile,status='new',form='formatted') write(6,1011)trim(gmtfile) 1011 format(6x,'mymedian5.f: Creating GMT batch file ',a) write(99,1030)trim(gmtfile) 1030 format('echo BEGIN batch file ',a) c ------------------------------------------------------------------ c - Open the input, pre-thinned, GMT-ready files c - VERTCON3: Coverage is coverage, whether for c - cd, v2 or rz data. Only going to plot "cd" coverage. c ------------------------------------------------------------------ c - Coverage (11-13): c gfncvacdlat = 'cvacdlat.'//trim(suffix1) c open(11,file=gfncvacdlat,status='old',form='formatted') c write(6,1010)trim(gfncvacdlat) c gfncvacdlon = 'cvacdlon.'//trim(suffix1) c open(12,file=gfncvacdlon,status='old',form='formatted') c write(6,1010)trim(gfncvacdlon) gfncvacdoht = 'cvacdoht.'//trim(suffix1) open(13,file=gfncvacdoht,status='old',form='formatted') write(6,1010)trim(gfncvacdoht) c - Vectors, Meters (21-24): c - Won't need to use the v2 or rz names unless 'v2special' is true, but going to c - establish them anyways. This is the philosophy of the fix c - on 2018 05 04 gfnvmacdoht = 'vmacdoht.'//trim(suffix1) gfnvmav2oht = 'vmav2oht.'//trim(suffix1) gfnvmarzoht = 'vmarzoht.'//trim(suffix1) c - If "v2special" we grid RZ data, not CD data, but we MUST carry the CD data c - along, thinning/dropping in parallel to what we do with the RZ data, since c - later (doit4 / checkgrid), we need access to the CD data to compare to the c - final grids. As such, open the "all" "CD" data here too, so we can perform c - identical tasks (thin/drop) on it that we perform on the "RZ" data... if(v2special)then open(23,file=gfnvmarzoht,status='old',form='formatted') write(6,1010)trim(gfnvmarzoht) open(92,file=gfnvmacdoht,status='old',form='formatted') write(6,1010)trim(gfnvmacdoht) else open(23,file=gfnvmacdoht,status='old',form='formatted') write(6,1010)trim(gfnvmacdoht) endif 1010 format(6x,'mymedian5.f: ', *'Opening existing unthinned file ',a) c ------------------------------------------------------------------ c - Open the output, thinned, GMT-ready files c - "t" meaning "thinned" c ------------------------------------------------------------------ c - Coverage (31-33): gfncvtcdoht = 'cvtcdoht.'//trim(suffix2) open(33,file=gfncvtcdoht,status='new',form='formatted') write(6,1022)trim(gfncvtcdoht) c - Vectors, Meters (41-44): gfnvmtcdoht = 'vmtcdoht.'//trim(suffix2) gfnvmtrzoht = 'vmtrzoht.'//trim(suffix2) gfnvmtv2oht = 'vmtv2oht.'//trim(suffix2) c - Updated on 2018/06/18: Although v2special still gets thinned on c - "RZ" data, we carry along a similarly thinned/dropped set c - of "CD" files too, for use in doit4.bat later (specifically checkgrid.f) c - If v2special, thin ON the RZ data, not CD data, BUT carry the thin CD data along in parallel if(v2special)then open(43,file=gfnvmtrzoht,status='new',form='formatted') write(6,1012)trim(gfnvmtrzoht) open(93,file=gfnvmtcdoht,status='new',form='formatted') write(6,1012)trim(gfnvmtcdoht) else open(43,file=gfnvmtcdoht,status='new',form='formatted') write(6,1012)trim(gfnvmtcdoht) endif 1012 format(6x,'mymedian5.f: ', *'Opening output thinned vector file ',a) 1022 format(6x,'mymedian5.f: ', *'Opening output thinned coverage file ',a) c ------------------------------------------------------------------ c - Open the output, dropped, GMT-ready files c - "d" meaning "dropped" c ------------------------------------------------------------------ c - Coverage (51-53): gfncvdcdoht = 'cvdcdoht.'//trim(suffix2) open(53,file=gfncvdcdoht,status='new',form='formatted') write(6,1023)trim(gfncvdcdoht) c - Vectors, Meters (61-64): gfnvmdcdoht = 'vmdcdoht.'//trim(suffix2) gfnvmdrzoht = 'vmdrzoht.'//trim(suffix2) gfnvmdv2oht = 'vmdv2oht.'//trim(suffix2) c - If v2special, thin ON the RZ data, not CD data, BUT carry the thin CD data along in parallel if(v2special)then open(63,file=gfnvmdrzoht,status='new',form='formatted') write(6,1013)trim(gfnvmdrzoht) open(94,file=gfnvmdcdoht,status='new',form='formatted') write(6,1013)trim(gfnvmdcdoht) else open(63,file=gfnvmdcdoht,status='new',form='formatted') write(6,1013)trim(gfnvmdcdoht) endif 1013 format(6x,'mymedian5.f: ', *'Opening output dropped vector file ',a) 1023 format(6x,'mymedian5.f: ', *'Opening output dropped coverage file ',a) c ------------------------------------------------------------------ c - Open the output, thinned, GMT-ready files for pushing c - through "surface" to get a grid. c - "t" meaning "thinned" c ------------------------------------------------------------------ c - Surface ready vectors, meters (71-73): gfnsmtcdoht = 'smtcdoht.'//trim(suffix2) gfnsmtrzoht = 'smtrzoht.'//trim(suffix2) gfnsmtv2oht = 'smtv2oht.'//trim(suffix2) c - If v2special, thin ON the RZ data, not CD data. However, as we are NOT going to c - run "SURFACE" on the thinned CD data (during v2special), we do NOT need to generate c - any parallel CD files for thin in "surface" style if(v2special)then open(73,file=gfnsmtrzoht,status='new',form='formatted') write(6,1014)trim(gfnsmtrzoht) else open(73,file=gfnsmtcdoht,status='new',form='formatted') write(6,1014)trim(gfnsmtcdoht) endif 1014 format(6x,'mymedian5.f: ', *'Opening output thinned vector file (for surface):',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,'mymedian5.f: Region= ',a,/, * 6x,'mymedian5.f: North = ',f12.6,/, * 6x,'mymedian5.f: South = ',f12.6,/, * 6x,'mymedian5.f: West = ',f12.6,/, * 6x,'mymedian5.f: East = ',f12.6) c ------------------------------------------------------- c - Get the header information necessary for a ".b" file c - This is the grid upon which thinning/dropping occurs. c - In "v2special" that is the RZ data (igridsec) c - In non-v2special that is the CD data but STILL is "igridsec" c - 1) Convert "igridsec" to decimal degrees c - 2) Count rows and columns c ------------------------------------------------------- dgla = dble(igridsec)/3600.d0 dglo = dble(igridsec)/3600.d0 nla=idnint((glamx-glamn)/dgla)+1 nlo=idnint((glomx-glomn)/dglo)+1 write(6,3001)glamn,glamx,glomn,glomx,dgla,dglo,nla,nlo 3001 format(6x,'mymedian5.f : Cell Structure:',/, *8x,'North = ',f16.10,/, *8x,'South = ',f16.10,/, *8x,'West = ',f16.10,/, *8x,'East = ',f16.10,/, *8x,'DLat = ',f16.10,/, *8x,'DLon = ',f16.10,/, *8x,'NLat = ',i16 ,/, *8x,'NLon = ',i16 ) c glamx=glamn+dgla*(nla-1) !*** register the far boundary c glomx=glomn+dglo*(nlo-1) c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c c - VERTCON 3 modification: No horizontal, so the first trip through c - the loop for "horizontal" is now gone. c c - Median filter in a loop. First loop is to process on absolute horizontal c - distance in METERS. Second loop on ellipsoid height, again on c - absolute horizontal value in METERS. (All "arcsecond" stuff c - is passively carried along to match the meter stuff. All the c - WORK is done in meters for median thinning). c - c - In the first loop, when done processing on horizontal distance, c - use the kept PIDs to generate the thinned LAT and thinned LON c - files (in this way, the PID list for both LAT and LON files will c - always be identical). Do NOT process medians on LAT and LON c - separately!!! c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c ----------------------------------------------------------------------- c nthinhor = 0 c ndrophor = 0 nthinoht = 0 ndropoht = 0 c do 2001 iloop = 1,2 do 2001 iloop = 2,2 if (iloop.eq.1)then lin = 24 ! vchor file c fused = gfnvmacdhor elseif(iloop.eq.2)then c - Here is the "all" file. In v2special it is "rz". Otherwise it is "cd". c - This is the data upon which medians and thin/drop decisions must be made. c - For v2special, these identical decisions must go in parallel to the c - "all/cd" file too, without actually analyzing the "all/cd" data itself. c - Just a perfect 1:1 correspondence. lin = 23 ! vcoht file if(v2special)then fused = gfnvmarzoht else fused = gfnvmacdoht endif endif write(6,2002)trim(fused) 2002 format(6x,'mymedian5.f: Median filtering file : ',a) 101 format(f16.10,1x,f15.10,1x,6x ,1x,12x ,1x,9x ,1x,f9.3,1x,a6) c - 2016 09 14: c2003 format(6x,'mymedian5.f: Point outside boundaries: ', c * f16.10,1x,f15.10,1x,f9.5,1x,a6) 2003 format(6x,'mymedian5.f: FATAL ERROR: ', * 'Point outside boundaries: ', * f16.10,1x,f15.10,1x,f9.5,1x,a6) c - Set the count of how many values are in each non-empty c - grid cell to zero. do 201 i=1,max ppcell(i)=0 201 continue ncells = 0 c ------------------------------------------------------------ c - Top of READ loop c - "ikt" is the count of all records in the vector file c ------------------------------------------------------------ ikt = 0 c write(6,*) ' before read loop' 100 read(lin,'(a)',end=777)card c - Note that we will work in (aka "the z value") METERS on both horizontal c - *AND* ellipsoid heights. read(card,101)glo,gla,z,pid ikt = ikt + 1 c - 2016 09 14: Changed this code...if there is ANY point that comes c - in from outside the grid, that will be a FATAL error c - in this program. if(gla.lt.glamn.or.gla.gt.glamx .or. * glo.lt.glomn.or.glo.gt.glomx) then write(6,2003)gla,glo,z,pid stop endif c - Store the read data cards(ikt)=card c - Store the values too zs(ikt) =z c - Set a default for ALL points that it will "not pass" c - (e.g. "be dropped") lpass(ikt)=.false. c - Determine which row and column our point falls in. Then c - convert that combination of row/col into a single c - value (encoded as ila*100000+ilo) ila=idnint((gla-glamn)/dgla)+1 ilo=idnint((glo-glomn)/dglo)+1 ipos=ila*100000+ilo c - Determine the distance our point is from the cell's center. c - This is used to break a tie when there are an even number c - of points in a cell. xla = glamn + (ila-1)*dgla xlo = glomn + (ilo-1)*dglo dist(ikt) = dsqrt((gla-xla)**2+(glo-xlo)**2) c - Is this point a new cell? do 202 j=1,ncells c - OLD cell c - Ratchet up the count of points in this cell c - Store the unique input number as being in this cell c - Jump out if(ipos.eq.poscode(j))then ppcell(j) = ppcell(j) + 1 whichcell(ikt) = j goto 203 endif 202 continue c - NEW cell (jump over this code from previous loop if c - an old cell were found) ncells = ncells + 1 poscode(ncells) = ipos ilapcell(ncells) = ila ilopcell(ncells) = ilo ppcell(ncells) = 1 c 2016 09 14: whichcell(ikt) = ncells c write(6,*) ' new cell : ',ncells,ipos 203 continue goto 100 777 continue c ------------------------------------------------------------ c - Bottom of READ loop c ------------------------------------------------------------ nkt = ikt write(6,2004) trim(fused),nkt,igridsec,ncells 2004 format( *6x,'mymedian5.f: Done Reading from : ',a,/, *6x,'mymedian5.f: Points in File : ',i10,/, *6x,'mymedian5.f: Cell Size (Arcseconds) : ',i10,/, *6x,'mymedian5.f: Number of Cells with Data : ',i10) write(6,2005) 2005 format(6x,'mymedian5.f: Begin median filtering') c ------------------------------------------------------------ c - TOP of SORTING section c ------------------------------------------------------------ c- Sort all of our data (nkt values) by their poscode write(6,2020) 2020 format(6x,'mymedian5.f: Sorting data...') c - 2016 09 14 write(6,744) nkt,max 744 format(6x,'mymedian5.f: Calling indexxi, nkt,max = ',i0,1x,i0) call indexxi(nkt,max,whichcell,indx) write(6,745) 745 format(6x,'mymedian5.f: Back from indexxi') c - 2016 09 14: Fixed, but kept commented out c - Spit out our data in WHICHCELL order c goto 400 c do 400 i=1,nkt c j = indx(i) c icell = whichcell(j) c if(mod(i,1000).eq.0)then c write(6,2010)cards(j),icell,ilapcell(icell), c * ilopcell(icell),poscode(icell) c endif c 400 continue 2010 format(i8,1x,a80,1x,4(i10)) c - Now go through each cell and sort it inumlo = 0 inumhi = 0 ikt = 0 do 401 icell=1,ncells inumlo = inumhi + 1 inumhi = inumlo + ppcell(icell) - 1 c write(6,*) ' inumlo/inumhi = ',inumlo,inumhi do 402 jpt=1,ppcell(icell) ikt = ikt + 1 j = indx(ikt) c write(6,2010)ikt,cards(j),icell,ilapcell(icell), c * ilopcell(icell),poscode(icell) 402 continue call getmedian(zs,max,inumlo,inumhi,indx, * maxppcell,dist,iiival) c - Set the median for this cell to "pass" (e.g. to be c - put in the "thinned" data file. if(iiival.lt.0)stop 10001 if(iiival.gt.nkt)stop 10002 lpass(iiival) = .true. 401 continue write(6,408) 408 format(6x,'mymedian5.f: Done finding medians',/, * 6x,'mymedian5.f: Populating thinned and dropped files') rewind(lin) ithin = 0 idrop = 0 c - 2016 09 14 do 451 i=1,nkt if(iloop.eq.1)then c read(11,'(a)')cardcvlat c read(12,'(a)')cardcvlon c read(21,'(a)')cardvmlat c read(22,'(a)')cardvmlon c read(24,'(a)')cardvmhor c read(26,'(a)')cardvslat c read(27,'(a)')cardvslon c read(29,'(a)')cardvshor c if(lpass(i))then c ithin = ithin + 1 c write(31,'(a)')cardcvlat c write(32,'(a)')cardcvlon c write(41,'(a)')cardvmlat c write(42,'(a)')cardvmlon c write(44,'(a)')cardvmhor c write(46,'(a)')cardvslat c write(47,'(a)')cardvslon c write(49,'(a)')cardvshor c - write the thinned vectors to the surface files c - NOTE: We need to extract LON/LAT/Value_in_correct_units c - from the overall "card", and the correct value c - depends on whether we're pulling out arcseconds or c - meters. c write(71,'(a,a,a)')cardvmlat(1:33),cardvmlat(64:72), c * cardvmlat(73:79) c write(72,'(a,a,a)')cardvmlon(1:33),cardvmlon(64:72), c * cardvmlat(73:79) c write(76,'(a,a,a)')cardvslat(1:33),cardvslat(54:62), c * cardvmlat(73:79) c write(77,'(a,a,a)')cardvslon(1:33),cardvslon(54:62), c * cardvmlat(73:79) c else c idrop = idrop + 1 c write(51,'(a)')cardcvlat c write(52,'(a)')cardcvlon c write(61,'(a)')cardvmlat c write(62,'(a)')cardvmlon c write(64,'(a)')cardvmhor c write(66,'(a)')cardvslat c write(67,'(a)')cardvslon c write(69,'(a)')cardvshor c endif else read(13,'(a)')cardcvoht c - File 23 is the "all" file, which was sorted. It is "cd" if not "v2special". c - It is "rz" if "v2special". Thus if v2special, we need to read (for carrying along) c - a "cd" file (unit 92) read(23,'(a)')cardvmoht if(v2special)read(92,'(a)')cardvmoht0 if(lpass(i))then ithin = ithin + 1 write(33,'(a)')cardcvoht write(43,'(a)')cardvmoht c - "Carry along" the thinned CD data under v2special: if(v2special)write(93,'(a)')cardvmoht0 c - write the thinned vectors to the surface file write(73,'(a,a,a)')cardvmoht(1:33),cardvmoht(64:72), * cardvmoht(73:79) else idrop = idrop + 1 write(53,'(a)')cardcvoht write(63,'(a)')cardvmoht c - "Carry along" the dropped CD data under v2special: if(v2special)write(94,'(a)')cardvmoht0 endif endif 451 continue write(6,746)ithin,idrop 746 format(6x,'mymedian5.f: Number Kept : ',i0,/, * 6x,'mymedian5.f: Number Dropped: ',i0) if(iloop.eq.1)then c nthinhor = ithin c ndrophor = idrop else nthinoht = ithin ndropoht = idrop endif c do 300 i=1,ncells c write(6,2006)i,poscode(i) c do 301 j=1,ppcell(i) c write(6,2007)j c 301 continue c 300 continue c ------------------------------------------------------------ c - BOTTOM of SORTING section c ------------------------------------------------------------ 2001 continue c ------------------------------------------------------- c - Based on everything we have so far, put a bunch c - of "surface" calls into the GMT batch file. c - These will be to turn the median-thinned data c - into grids. DO NOT GENERATE A CALL TO SURFACE c - IF THE NUMBER OF THINNED DATA IS ZERO. c - Because "surface", as part of GMT, returns its c - own binary grid format (called ".grd" herein), it must be c - converted to ".b" format. The easiest way to c - do so is to use the GMT built-in routine "grd2xyz" c - with the "-bos" extension to create a binary XYZ c - list file. Then run Dennis's homemade "xyz2b.for" c - routine to finally arrive at the ".b" binary grid format. c c - 02018 03 28: GMT5: -bos is replaced with -bof c c ------------------------------------------------------- cmidlat=dcos(((glamn+glamx)/2.d0)*d2r) rfnvmtcdoht00 = 'vmtcdoht.'//trim(suffix2t00)//'.grd' rfnvmtcdoht04 = 'vmtcdoht.'//trim(suffix2t04)//'.grd' rfnvmtcdoht10 = 'vmtcdoht.'//trim(suffix2t10)//'.grd' rfnvmtrzoht00 = 'vmtrzoht.'//trim(suffix2t00)//'.grd' rfnvmtrzoht04 = 'vmtrzoht.'//trim(suffix2t04)//'.grd' rfnvmtrzoht10 = 'vmtrzoht.'//trim(suffix2t10)//'.grd' c if(v2special)then c else c endif sadrfnvmtcdoht1000 = 'vmtcdoht.'//trim(suffix2d3)//'.grd' sadrfnvmtrzoht1000 = 'vmtrzoht.'//trim(suffix2d3)//'.grd' c if(v2special)then c else c endif zfnvmtcdoht00 = 'vmtcdoht.'//trim(suffix2t00)//'.xyz' zfnvmtcdoht04 = 'vmtcdoht.'//trim(suffix2t04)//'.xyz' zfnvmtcdoht10 = 'vmtcdoht.'//trim(suffix2t10)//'.xyz' zfnvmtrzoht00 = 'vmtrzoht.'//trim(suffix2t00)//'.xyz' zfnvmtrzoht04 = 'vmtrzoht.'//trim(suffix2t04)//'.xyz' zfnvmtrzoht10 = 'vmtrzoht.'//trim(suffix2t10)//'.xyz' c if(v2special)then c else c endif bfnvmtcdoht00 = 'vmtcdoht.'//trim(suffix2t00)//'.b' bfnvmtcdoht04 = 'vmtcdoht.'//trim(suffix2t04)//'.b' bfnvmtcdoht10 = 'vmtcdoht.'//trim(suffix2t10)//'.b' bfnvmtrzoht00 = 'vmtrzoht.'//trim(suffix2t00)//'.b' bfnvmtrzoht04 = 'vmtrzoht.'//trim(suffix2t04)//'.b' bfnvmtrzoht10 = 'vmtrzoht.'//trim(suffix2t10)//'.b' c if(v2special)then c else c endif c - Difference grid names dbfnvmtcdoht1000 = 'vmtcdoht.'//trim(suffix2d1)//'.b' adbfnvmtcdoht1000 = 'vmtcdoht.'//trim(suffix2d2)//'.b' sadbfnvmtcdoht1000 = 'vmtcdoht.'//trim(suffix2d3)//'.b' dbfnvmtrzoht1000 = 'vmtrzoht.'//trim(suffix2d1)//'.b' adbfnvmtrzoht1000 = 'vmtrzoht.'//trim(suffix2d2)//'.b' sadbfnvmtrzoht1000 = 'vmtrzoht.'//trim(suffix2d3)//'.b' c if(v2special)then c else c endif c - Due to a bug in GMT, the call will use "-I*m" rather c - than "-I*s". As such, convert igridsec to xgridmin. xgridmin = dble(igridsec)/60.d0 c - 2018 05 16: Two grid spacings: xgridmin2 = dble(igridsec2)/60.d0 c - Do NOT generate a call to surface if there is no data c - to grid. c - Call to grid the thinned latitudes if(nthinhor.ne.0)then c write(99,501)trim(gfnsmtcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlat00),0.0,cmidlat c write(99,501)trim(gfnsstcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlat00),0.0,cmidlat c write(99,501)trim(gfnsmtcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlat04),0.4,cmidlat c write(99,501)trim(gfnsstcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlat04),0.4,cmidlat c write(99,501)trim(gfnsmtcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlat10),1.0,cmidlat c write(99,501)trim(gfnsstcdlat),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlat10),1.0,cmidlat endif c - Call to grid the thinned longitudes if(nthinhor.ne.0)then c write(99,501)trim(gfnsmtcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlon00),0.0,cmidlat c write(99,501)trim(gfnsstcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlon00),0.0,cmidlat c write(99,501)trim(gfnsmtcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlon04),0.4,cmidlat c write(99,501)trim(gfnsstcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlon04),0.4,cmidlat c write(99,501)trim(gfnsmtcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvmtcdlon10),1.0,cmidlat c write(99,501)trim(gfnsstcdlon),glomn,glomx,glamn,glamx, c * xgridmin,trim(rfnvstcdlon10),1.0,cmidlat endif c - Call to grid the thinned orthometric heights c - Note, these are RZ under V2special or CD under others if(nthinoht.ne.0)then if(v2special)then write(99,501)trim(gfnsmtrzoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtrzoht00),0.0,cmidlat write(99,501)trim(gfnsmtrzoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtrzoht04),0.4,cmidlat write(99,501)trim(gfnsmtrzoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtrzoht10),1.0,cmidlat else write(99,501)trim(gfnsmtcdoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtcdoht00),0.0,cmidlat write(99,501)trim(gfnsmtcdoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtcdoht04),0.4,cmidlat write(99,501)trim(gfnsmtcdoht),glomn,glomx,glamn,glamx, * xgridmin,trim(rfnvmtcdoht10),1.0,cmidlat endif endif c - Calls to convert ".grd" to ".xyz" if(nthinhor.ne.0)then c write(99,502)trim(rfnvmtcdlat00),trim(zfnvmtcdlat00) c write(99,502)trim(rfnvstcdlat00),trim(zfnvstcdlat00) c write(99,502)trim(rfnvmtcdlat04),trim(zfnvmtcdlat04) c write(99,502)trim(rfnvstcdlat04),trim(zfnvstcdlat04) c write(99,502)trim(rfnvmtcdlat10),trim(zfnvmtcdlat10) c write(99,502)trim(rfnvstcdlat10),trim(zfnvstcdlat10) endif if(nthinhor.ne.0)then c write(99,502)trim(rfnvmtcdlon00),trim(zfnvmtcdlon00) c write(99,502)trim(rfnvstcdlon00),trim(zfnvstcdlon00) c write(99,502)trim(rfnvmtcdlon04),trim(zfnvmtcdlon04) c write(99,502)trim(rfnvstcdlon04),trim(zfnvstcdlon04) c write(99,502)trim(rfnvmtcdlon10),trim(zfnvmtcdlon10) c write(99,502)trim(rfnvstcdlon10),trim(zfnvstcdlon10) endif if(nthinoht.ne.0)then if(v2special)then write(99,502)trim(rfnvmtrzoht00),trim(zfnvmtrzoht00) write(99,502)trim(rfnvmtrzoht04),trim(zfnvmtrzoht04) write(99,502)trim(rfnvmtrzoht10),trim(zfnvmtrzoht10) else write(99,502)trim(rfnvmtcdoht00),trim(zfnvmtcdoht00) write(99,502)trim(rfnvmtcdoht04),trim(zfnvmtcdoht04) write(99,502)trim(rfnvmtcdoht10),trim(zfnvmtcdoht10) endif endif c - Calls to convert ".xyz" to ".b" if(nthinhor.ne.0)then c write(99,503)trim(zfnvmtcdlat00),trim(bfnvmtcdlat00) c write(99,503)trim(zfnvstcdlat00),trim(bfnvstcdlat00) c write(99,503)trim(zfnvmtcdlat04),trim(bfnvmtcdlat04) c write(99,503)trim(zfnvstcdlat04),trim(bfnvstcdlat04) c write(99,503)trim(zfnvmtcdlat10),trim(bfnvmtcdlat10) c write(99,503)trim(zfnvstcdlat10),trim(bfnvstcdlat10) endif if(nthinhor.ne.0)then c write(99,503)trim(zfnvmtcdlon00),trim(bfnvmtcdlon00) c write(99,503)trim(zfnvstcdlon00),trim(bfnvstcdlon00) c write(99,503)trim(zfnvmtcdlon04),trim(bfnvmtcdlon04) c write(99,503)trim(zfnvstcdlon04),trim(bfnvstcdlon04) c write(99,503)trim(zfnvmtcdlon10),trim(bfnvmtcdlon10) c write(99,503)trim(zfnvstcdlon10),trim(bfnvstcdlon10) endif if(nthinoht.ne.0)then if(v2special)then write(99,503)trim(zfnvmtrzoht00),trim(bfnvmtrzoht00) write(99,503)trim(zfnvmtrzoht04),trim(bfnvmtrzoht04) write(99,503)trim(zfnvmtrzoht10),trim(bfnvmtrzoht10) else write(99,503)trim(zfnvmtcdoht00),trim(bfnvmtcdoht00) write(99,503)trim(zfnvmtcdoht04),trim(bfnvmtcdoht04) write(99,503)trim(zfnvmtcdoht10),trim(bfnvmtcdoht10) endif endif 501 format( *'surface ',a,' -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s, *' -I',f0.2,'m -G',a,' -T',f3.1,' -A',s,f6.4,' -C0.01 -V') c 501 format( c *'surface ',a,' -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s, c *' -I',f0.2,'m -G',a,' -T0.4 -A',s,f6.4,' -C0.01 -V') c - GMT5 change: c 502 format( c *'grd2xyz ',a,' -bos > ',a) 502 format( *'grd2xyz ',a,' -bof > ',a) 503 format( *'xyz2b << !',/, *a,/,a,/,'!') 504 format( *'subtrc << !',/, *a,/,a,/,a,/,'!') 505 format( *'gabs << !',/, *a,/,a,/,'!') 506 format( *'gscale << !',/, *a,/,f10.5,/,a,/,'!') c - 2018 03 28: GMT5: -bis becomes -bif 507 format( *'b2xyz << !',/,a,/,'!',/, *'xyz2grd temp.xyz -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s, *' -I',f0.2,'m -bif -G',a,/, *'rm -f temp.xyz') c 507 format( c *'b2xyz << !',/,a,/,'!',/, c *'xyz2grd temp.xyz -R',f9.5,'/',f9.5,'/',sp,f9.5,'/',f9.5,s, c *' -I',f0.2,'m -bis -G',a,/, c *'rm -f temp.xyz') c - New, as of 10/27/2015: c - Subtract the T=0.0 transformation grid c - from the T=1.0 transformation grid. c - Then take the absolute value of that c - difference, and finally scale that c - grid by 1/1.6929. Do this for c - all 5 possible transformation grids c - All of these calls should be in the c - "gmtbat02" file. c - Differences first: if(nthinhor.ne.0)then c write(99,504)trim(bfnvmtcdlat10),trim(bfnvmtcdlat00), c * trim(dbfnvmtcdlat1000) c write(99,504)trim(bfnvstcdlat10),trim(bfnvstcdlat00), c * trim(dbfnvstcdlat1000) c write(99,504)trim(bfnvmtcdlon10),trim(bfnvmtcdlon00), c * trim(dbfnvmtcdlon1000) c write(99,504)trim(bfnvstcdlon10),trim(bfnvstcdlon00), c * trim(dbfnvstcdlon1000) endif if(nthinoht.ne.0)then if(v2special)then write(99,504)trim(bfnvmtrzoht10),trim(bfnvmtrzoht00), * trim(dbfnvmtrzoht1000) else write(99,504)trim(bfnvmtcdoht10),trim(bfnvmtcdoht00), * trim(dbfnvmtcdoht1000) endif endif c - Make them absolute values: if(nthinhor.ne.0)then c write(99,505)trim(dbfnvmtcdlat1000),trim(adbfnvmtcdlat1000) c write(99,505)trim(dbfnvstcdlat1000),trim(adbfnvstcdlat1000) c write(99,505)trim(dbfnvmtcdlon1000),trim(adbfnvmtcdlon1000) c write(99,505)trim(dbfnvstcdlon1000),trim(adbfnvstcdlon1000) endif if(nthinoht.ne.0)then if(v2special)then write(99,505)trim(dbfnvmtrzoht1000),trim(adbfnvmtrzoht1000) else write(99,505)trim(dbfnvmtcdoht1000),trim(adbfnvmtcdoht1000) endif endif c - Finally, scale by 1/1.6929 if(nthinhor.ne.0)then c write(99,506)trim(adbfnvmtcdlat1000),rng2std, c * trim(sadbfnvmtcdlat1000) c write(99,506)trim(adbfnvstcdlat1000),rng2std, c * trim(sadbfnvstcdlat1000) c write(99,506)trim(adbfnvmtcdlon1000),rng2std, c * trim(sadbfnvmtcdlon1000) c write(99,506)trim(adbfnvstcdlon1000),rng2std, c * trim(sadbfnvstcdlon1000) endif if(nthinoht.ne.0)then if(v2special)then write(99,506)trim(adbfnvmtrzoht1000),rng2std, * trim(sadbfnvmtrzoht1000) else write(99,506)trim(adbfnvmtcdoht1000),rng2std, * trim(sadbfnvmtcdoht1000) endif endif c - After all this, I now have the "d3" files as ".b" grids, c - but in order to plot them, I'll have to convert them c - to ".grd" files. So here we go... if(nthinhor.ne.0)then c write(99,507)trim(sadbfnvmtcdlat1000) c * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvmtcdlat1000) c write(99,507)trim(sadbfnvstcdlat1000) c * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvstcdlat1000) c write(99,507)trim(sadbfnvmtcdlon1000) c * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvmtcdlon1000) c write(99,507)trim(sadbfnvstcdlon1000) c * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvstcdlon1000) endif c - Converting the "d3" .b grid to .grd: c - In "v2special", the d3 grid is RZ data c - Otherwise the d3 grid is CD data if(nthinoht.ne.0)then if(v2special)then write(99,507)trim(sadbfnvmtrzoht1000) * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvmtrzoht1000) else write(99,507)trim(sadbfnvmtcdoht1000) * ,glomn,glomx,glamn,glamx,xgridmin,trim(sadrfnvmtcdoht1000) endif endif c --------------------------------------------------------------------------------- c --------------------------------------------------------------------------------- c - VERTCON 3: NGVD29/NAVD88/CONUS ("v2special") issue (BEGIN): c --------------------------------------------------------------------------------- c --------------------------------------------------------------------------------- c - Updated 2018 05 16 to allow for TWO grid spacings. One for thinning/gridding c - the RZ data, and one for the grid spacing of the final CD grid. To achieve c - CD grids on this 2nd spacing, need to: c 1) Run "regrd2" on the RZ gridded data (at "igridsec") to get it on "igridsec2" c 2) Find the V2 grid at "igridsec2" c 3) Add #1 and #2 together to get CD data on "igridsec2" c - Finally here at the end, if I've been thinning/dropping/gridding c - on the "rz" files (because I'm under "v2special", which is c - the Remove/Compute/Restore process set up specially built c - around the VERTCON 2 grid), it's time to add back the c - "v2" grid itself, to yield the following: c a) An actual "cd" grid in ".b" format c b) The "cd" grid in ".grd" format for plotting. c - See DRU-12, p. 129 for details. c - See DRU-12, p. 133-4 for more about 2 grid spacings if(v2special)then c - First, get the gridded "RZ" data from "igridsec" spacing c - to "igridsec2" spacing... c - The "rz" transformation grid file in .b form at "igridsec" = bfnvmtrzoht04 c - The "rz" transformation grid file in .b form at "igridsec2" = bfnvmtrzoht04x = "temp.b" (to be killed) c - As we will kill this file, let's just give it a temp name bfnvmtrzoht04x = 'temp.b' c - Convert old rows/cols counts to new dglax = dble(igridsec2)/3600.d0 dglox = dble(igridsec2)/3600.d0 nlax=idnint((glamx-glamn)/dglax)+1 nlox=idnint((glomx-glomn)/dglox)+1 c - Call "regrd2" to do the re-gridding to "igridsec2" grid levels write(99,751)trim(bfnvmtrzoht04),trim(bfnvmtrzoht04x), * nlax,nlox 751 format('regrd2 << ! ',/, * a,/,a,/,i0,/,i0,/,'!') c - Determine the right "vertcon2" grid to use for this grid spacing: write(dumgridsec2,'(i4.4)')igridsec2 c v2fnameb = '/home/dru/VERTCON/VERTCON2/' c * //'vertcon2.0.v3bounds.sec'//dumgridsec2// c * '.meters.b' c v2fnameg = '/home/dru/VERTCON/VERTCON2/' c * //'vertcon2.0.v3bounds.sec'//dumgridsec2// c * '.meters.grd' c - Changed to a directory BELOW the /home/dru/VERTCON3 level c - as per conversation with ALB on 12/20/2018, to assist c - in simplifying the creation of the digital archive v2fnameb = 'VERTCON2/' * //'vertcon2.0.v3bounds.sec'//dumgridsec2// * '.meters.b' v2fnameg = 'VERTCON2/' * //'vertcon2.0.v3bounds.sec'//dumgridsec2// * '.meters.grd' c write(dumgridsec,'(i4.4)')igridsec c v2fnameb = '/home/dru/VERTCON/VERTCON2/' c * //'vertcon2.0.v3bounds.sec'//dumgridsec// c * '.meters.b' c v2fnameg = '/home/dru/VERTCON/VERTCON2/' c * //'vertcon2.0.v3bounds.sec'//dumgridsec// c * '.meters.grd' c - Now call an "addem" to add the above "temp.b" file c - to the appropriate, pre-existing V2 grid at "igridsec2" c - (Note that these pre-existing grids were created outside of this c - batch file process, and are all derived from the original c - 3' grid released as "VERTCON 2". See DRU-12, p. 129) c c - Add the thinned/gridded/T=0.4/rz ".b" file, which has now been c - regridded to "igridsec2" levels (and is now called "temp.b") to the c - proper "v2" grid: c - Add the "rz" transformation grid file in .b form at "igridsec2" (bfnvmtrzoht04x) c - to the "v2" transformation grid file in .b form at "igridsec2" (v2fnameb) c - to get the "cd" transformation grid file in .b form at "igridsec2" (bfnvmtcdoht04) c - "addem" first: write(99,803) * trim(bfnvmtrzoht04x), * trim(v2fnameb), * trim(bfnvmtcdoht04), * trim(bfnvmtrzoht04x), * trim(v2fnameb), * trim(bfnvmtcdoht04), * trim(bfnvmtrzoht04x), * trim(v2fnameb), * trim(bfnvmtcdoht04) c write(99,803) c * trim(bfnvmtrzoht04), c * trim(v2fnameb), c * trim(bfnvmtcdoht04), c * trim(bfnvmtrzoht04), c * trim(v2fnameb), c * trim(bfnvmtcdoht04), c * trim(bfnvmtrzoht04), c * trim(v2fnameb), c * trim(bfnvmtcdoht04) c - Then "b2xyz" and "xyz2grd" the "cd" file of "04" style: c - To get it into ".grd" form, needed for color plots with GMT write(99,507)trim(bfnvmtcdoht04) * ,glomn,glomx,glamn,glamx,xgridmin2,trim(rfnvmtcdoht04) c write(99,507)trim(bfnvmtcdoht04) c * ,glomn,glomx,glamn,glamx,xgridmin,trim(rfnvmtcdoht04) c- Aaaaaaaand, additionally, I need to create a dummy version of the existing VERTCON 2 grid c - See DRU-12, p. 130. write(99,509)trim(v2fnameg),'vmtv2oht.ngvd29.navd88.conus.'// * filter//'.'//trim(agridsec)//'.'//trim(agridsec2)//'.04.grd' c write(99,509)trim(v2fnameg),'vmtv2oht.ngvd29.navd88.conus.'// c * filter//'.'//trim(agridsec)//'.04.grd' write(99,509)trim(v2fnameb),'vmtv2oht.ngvd29.navd88.conus.'// * filter//'.'//trim(agridsec)//'.'//trim(agridsec2)//'.04.b' c write(99,509)trim(v2fnameb),'vmtv2oht.ngvd29.navd88.conus.'// c * filter//'.'//trim(agridsec)//'.04.b' 509 format('cp ',a,' ',a) c - Finally, kill the "temp.b" grid write(99,510)trim(bfnvmtrzoht04x) 510 format('rm -f ',a) endif c --------------------------------------------------------------------------------- c --------------------------------------------------------------------------------- c - VERTCON 3: NGVD29/NAVD88/CONUS ("v2special") issue (END): c --------------------------------------------------------------------------------- c --------------------------------------------------------------------------------- 803 format( *'# ------------------------------',/, *'# Adding ',a,' and ',a,' = ',a,/, *'# ------------------------------',/, *'echo Adding ',a,' and ',a,' = ',a,/, *'addem << !',/, *a,/,a,/,a,/,'!') write(99,1031)trim(gmtfile) 1031 format('echo END batch file ',a) close(99) write(6,9999) 9999 format('END mymedian5.f') 2006 format('Cell ',i8,' Poscode: ',i0) 2007 format(6x,'Pt : ',i5) end c c ----------------------------------------------- c subroutine getmedian(zs,max,inumlo,inumhi,indx, *maxppcell,dist,iiival) real*8 zs(max),dist(max) integer*4 indx(max),indx2(maxppcell) integer*4 inumlo,inumhi real*8 zz(maxppcell) c write(6,1)inumlo,inumhi c 1 format('Inside getmedian',/, c *'inumlo,inumhi = ',i8,1x,i8) c - When the whole set of "zs" data is sorted by c - "whichcell", then the "indx" values are our c - key to that sort. The "zs" data remains c - in "read order" (1-nkt), but the indx values c - (also in "read order" point to the sorted c - order. c write(6,*) ' UNSORTED data: ' nval = inumhi-inumlo+1 iq = 0 do 2 i=inumlo,inumhi iq = iq + 1 zz(iq) = zs(indx(i)) c write(6,3)i,iq,zz(iq) 2 continue 3 format(i10,1x,i10,1x,f20.10) c - Sort this mini vector c call hpsort(nval,maxppcell,zz) call indexxd(nval,maxppcell,zz,indx2) c write(6,*) ' SORTED data: ' c do 4 i=1,nval c write(6,3)i,indx2(i),zz(indx2(i)) c 4 continue 5 format(i10,1x,f20.10) c - True median comes from an odd number of values. If c - it is even, take the one which sits closest to the c - center of the cell. if(mod(nval,2).eq.1)then ival = (nval+1)/2 else ival1 = nval/2 ival2 = ival1 + 1 c write(6,*) ' MedianLo = ',zz(indx2(ival1)) c write(6,*) ' MedianHi = ',zz(indx2(ival2)) c write(6,*) ' Dist Lo = ',dist(indx2(ival1)) c write(6,*) ' Dist Hi = ',dist(indx2(ival2)) if(dist(indx2(ival1)).lt.dist(indx2(ival2)))then c write(6,*) ' Median Set to LO = ',zz(indx2(ival1)) ival = ival1 else c write(6,*) ' Median Set to HI = ',zz(indx2(ival2)) ival = ival2 endif endif c - Note that "iiival" is the read-order number of the c - point chosen as median for the cell we're analyzing c - in this subroutine. c write(6,*) ' Median = ',zz(indx2(ival)) igood1 = indx2(ival) c write(6,*) ' ival = ',ival c write(6,*) ' indx2(ival) = ',indx2(ival) iival = inumlo + indx2(ival) - 1 c write(6,*) ' iival = ',iival iiival = indx(iival) c write(6,*) ' indx(iival) = ',indx(iival) return end include 'Subs/getgridbounds.f' include '/home/dru/NumRec/Modified/indexxi.for' include '/home/dru/NumRec/Double/indexxd.for'