c - Program intg, Version 1.3, August 29, 2003 c - "Interpolate from Geoid files" c - Program will do the following: c - 1) Take one or more associated geoid files c - in binary (*.bin) format and compute, via c - biquadratic interpolation, the c - geoid undulation at any given c - latitude/longitude combination. c - This interpolation can be done for c - one or more points, interactively c - or through input/output files. c - Supported Geoid models as of version 1.3: c 1) GEOID03 3) GEOID99 c 2) USGG2003 4) G99SSS c - Note that "*.bin" format is binary, unformatted, direct access c - and that the order of bytes depends on which platform the c - file was created. c c VERSION DATE PRIMARY CONTACT c 1.0 Sep 22, 1999 Dru A. Smith c 1.1 Oct 7, 1999 Dru A. Smith c 1.2 Nov 2, 1999 Dru A. Smith c 1.3 Aug 29, 2003 Dan Roman c - Modifications for version 1.1: c 1) Fixed a bug that caused Blue Book output c files to fill up with repeated data c unendingly. c 2) Added code to allow users to enter longitude c in positive west coordinates if they want. c 3) Fixed a bug to allow negative signs to c be used in latitude & longitude entries c 4) Fixed a bug which cut off the left c byte (This is a DOS/Unix issue) c 5) Added more descriptive text about c allowable formats for input values c 6) Fixed a bug that didn't allow 81,82 nor 85 records c of the Blue Book to pass through c 7) Fixed a bug that computed "xx" and "yy" wrong in c the interpolation routine. c c - Modifications for version 1.2: c 1) Fixed the double-backslash typo that appears in c the "PC example" directory. c 2) Added a way to exit the program after format descriptions c are given. c c - Modifications for version 1.3: c 1) Added in extar options for models (usgg2003, c ushg2003, and xxusg). c 2) Started work on checking for platform configuration c and byte swap. c note) bug fix 2004apr16, at lines 758, 834, and 904/2736 c see Soft_Req #1726 c c c For further information, questions, or comments: c Daniel R. Roman, Ph.D. c NOAA, National Geodetic Survey, N/NGS6 c U.S.A. c Phone : 301-713-3202 c Fax : 301-713-4172 c e-mail : dan.roman@noaa.gov *********************************************************************** * * * DISCLAIMER * * * * THIS PROGRAM AND SUPPORTING INFORMATION IS FURNISHED BY THE * * GOVERNMENT OF THE UNITED STATES OF AMERICA, AND IS ACCEPTED AND * * USED BY THE RECIPIENT WITH THE UNDERSTANDING THAT THE UNITED STATES * * GOVERNMENT MAKES NO WARRANTIES, EXPRESS OR IMPLIED, CONCERNING THE * * ACCURACY, COMPLETENESS, RELIABILITY, OR SUITABILITY OF THIS * * PROGRAM, OF ITS CONSTITUENT PARTS, OR OF ANY SUPPORTING DATA. * * * * THE GOVERNMENT OF THE UNITED STATES OF AMERICA SHALL BE UNDER NO * * LIABILITY WHATSOEVER RESULTING FROM ANY USE OF THIS PROGRAM. THIS * * PROGRAM SHOULD NOT BE RELIED UPON AS THE SOLE BASIS FOR SOLVING A * * PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO PERSON * * OR PROPERTY. * * * * THIS PROGRAM IS PROPERTY OF THE GOVERNMENT OF THE UNITED STATES * * OF AMERICA. THEREFORE, THE RECIPIENT FURTHER AGREES NOT TO ASSERT * * PROPRIETARY RIGHTS THEREIN AND NOT TO REPRESENT THIS PROGRAM TO * * ANYONE AS BEING OTHER THAN A GOVERNMENT PROGRAM. * * * *********************************************************************** character*256 dirnam character*256 ofnam,ifnam character*1 ofyn,dumyn logical outfil logical poseast character*13 pcdirex integer*4 lin(50) character*256 fnam(50) integer*4 ios(50) real*8 glamn(50),glamx(50),glomn(50),glomx(50) real*8 dla(50),dlo(50) integer*4 nla(50),nlo(50) integer*4 ikind(50) real*4 val character*80 card,card2,card86 character*40 b40,sname character*40 xlatc,xlonc real*8 xlat,xlon logical ne,se,we,ee c variables and functions for handling endian problem logical reverse_bytes(50) common /endian/reverse_bytes real*8 frev8 c real*4 frev4 integer*4 irev4 c - Variables for statistics of the run integer*4 kount,keep real*8 max,min,maxlat,maxlon,minlat,minlon real*8 xn,xs,xw,xe real*8 ave,std,rms c - Variables to get around the recl=4 issue for our real*8 c - header variables real*4 glamnx(2),glomnx(2),dlax(2),dlox(2) real*8 glamny,glomny,dlay,dloy equivalence(glamnx(1),glamny) equivalence(glomnx(1),glomny) equivalence(dlax(1),dlay) equivalence(dlox(1),dloy) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Use common statements to cut down on the c - amount of data sent through 'call' statements ccccccccccccccccccccccccccccccccccccccccccccccccccc common /grid1/ glamn,glamx,glomn,glomx common /grid2/ dla,dlo,nla,nlo,ikind common /grid3/ ios common /edges/ ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Zero out the statistics for this run ccccccccccccccccccccccccccccccccccccccccccccccccccc ave = 0.d0 std = 0.d0 rms = 0.d0 keep = 0 kount = 0 max = -99999999.d0 min = +99999999.d0 xn = -99999.d0 xs = +99999.d0 xw = 99999.d0 xe = -99999.d0 b40 = ' '// * ' ' ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Write out the introductory/disclaimer screens ccccccccccccccccccccccccccccccccccccccccccccccccccc call intro ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Define the file numbers c - linn = Input file of lat/lon values c - lout = Output file of lat/lon/geoid values c - lin = set of geoid files (*.bin) ccccccccccccccccccccccccccccccccccccccccccccccccccc do 1 i=1,50 lin(i) = 10 + i 1 continue linn = 1 lout = 10 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Which geoid? ccccccccccccccccccccccccccccccccccccccccccccccccccc 108 write(6,101) 101 format(/,1x,70('-'),/, * ' Which geoid model do you wish to use?',//, * ' 1 = GEOID03 2 = USGG2003',/, * ' 3 = GEOID99 4 = G99SSS',/, * ' 5 = Experimental Geoid (XUSHG) ',/, * ' ',/, * ' 99 = END PROGRAM',//, * ' -> ',$) read(5,*)imodel if(imodel.eq.99)stop if(imodel.lt.1 .or. imodel.gt.5)goto 108 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Which directory are the geoid files in? ccccccccccccccccccccccccccccccccccccccccccccccccccc c - For the Unix Library Copy, I set the directory c location for the user. Modification made by c Bessie Thompson. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(imodel .eq. 1) then dirnam = '/ngslib/data/Geoid/Geoid03/' elseif(imodel .eq. 2) then dirnam = '/ngslib/data/Geoid/USGG2003/' elseif(imodel.eq.3) then dirnam = '/ngslib/data/Geoid/Geoid99/' elseif(imodel.eq.4) then dirnam = '/ngslib/data/Geoid/G99sss/' elseif(imodel.eq.5) then c dirnam = '/ngslib/data/Geoid/Geoid99/' endif c c pcdirex = char(39)//'C:'//char(92)//'GEOID03'//char(92)//char(39) c write(6,102)pcdirex c 102 format(/,1x,70('-'),/, c * ' What is the **FULL** directory name (including', c * ' trailing slashes) ',/, c * ' where the geoid (*.bin) files may be', c * ' found?',/, c * 5x,' (Unix Example: ''/export/home/geoid99/'') ',/, c * 5x,' (PC Example: ',a13',) ',/, c * 5x,' Hit to default to the current directory',//, c * ' -> ',$) c read(5,'(a)')dirnam ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Create the list of files that must be opened, c - and open them. Return with a count of how c - many files were opened, and a flag (ios) c - of which files are open and which are not. ccccccccccccccccccccccccccccccccccccccccccccccccccc call files(imodel,dirnam,nfiles,fnam,lin,nff) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Read the headers of all geoid files which c - where opened, and store that information. c - Compute the max lat/lon from the header information. ccccccccccccccccccccccccccccccccccccccccccccccccccc do 115 i=1,nfiles if(ios(i).eq.0)then read(lin(i),rec= 1)glamnx(1) read(lin(i),rec= 2)glamnx(2) read(lin(i),rec= 3)glomnx(1) read(lin(i),rec= 4)glomnx(2) read(lin(i),rec= 5)dlax(1) read(lin(i),rec= 6)dlax(2) read(lin(i),rec= 7)dlox(1) read(lin(i),rec= 8)dlox(2) read(lin(i),rec= 9)nla(i) read(lin(i),rec=10)nlo(i) read(lin(i),rec=11)ikind(i) C glamn(i) = glamny glomn(i) = glomny dla(i) = dlay dlo(i) = dloy C c convert big endian values and to little endian if(ikind(i).eq.1) then reverse_bytes(i)=.FALSE. elseif(ikind(i).eq.irev4(1)) then reverse_bytes(i)=.TRUE. else write(6,6115) fnam(i) 6115 format('The file ', A, 'is in an unrecognized format.'/ * 'intg is terminating') stop endif if(reverse_bytes(i)) then glamn(i)= frev8(glamn(i)) glomn(i)= frev8 (glomn(i)) dla(i) = frev8 (dla(i)) dlo(i) = frev8 (dlo(i)) nla(i) = irev4 (nla(i)) nlo(i) = irev4 (nlo(i)) ikind(i)=irev4 (ikind(i)) endif c glamx(i) = glamn(i) + (nla(i)-1)*dla(i) glomx(i) = glomn(i) + (nlo(i)-1)*dlo(i) endif 115 continue ccccccccccccccccccccccccccccccccccccccccccccccccccc c - How to input? ccccccccccccccccccccccccccccccccccccccccccccccccccc 110 write(6,103) 103 format(/,1x,70('-'),/, * ' How would you like to input the data? ',/, * ' 1 = By Keyboard (with prompts)',/, * ' 2 = By File (using allowed formats)',//, * ' -> ',$) read(5,*)iinput if(iinput.lt.1 .or. iinput.gt.2)goto 110 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If using an input file, find out which format c - it is in. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iinput.eq.2)then 105 write(6,104) 104 format(/,1x,70('-'),/, * ' Which format will you use for input? ',/, * ' 1 = Free Format (For Geoid) Type 1',/, * ' 2 = Free Format (For Geoid) Type 2',/, * ' 3 = NGS Blue Book Format (*80* and *86* records)',//, * ' 0 = End Program',/, * ' 99 = Please explain to me the formats',//, * ' -> ',$) read(5,*)iform if(iform.eq.0)stop ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Explain the formats if it was requested ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iform.eq.99)then call expform goto 105 endif endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Which positive longitude convention? c - DON'T ask this if the input is Blue Book ccccccccccccccccccccccccccccccccccccccccccccccccccc if(.not.(iinput.eq.2 .and. iform.eq.3))then 1103 write(6,1102) 1102 format(/,1x,70('-'),/, * ' Which longitude convention will you use? ',/, * ' 1 - Positive EAST',/, * ' 2 - Positive WEST',//, * ' -> ',$) read(5,*)ipos if(ipos.lt.1 .or. ipos.gt.2)goto 1103 if(ipos.eq.1)poseast = .true. if(ipos.eq.2)poseast = .false. endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Open the Input file if necessary ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iinput.eq.2)then 112 write(6,109) 109 format(/,1x,70('-'),/, * ' What is your input file name? ',//, * ' -> ',$) read(5,'(a)')ifnam open(linn,file=ifnam,status='old',form='formatted', * iostat=ioss) if(ioss.ne.0)then write(6,111) 111 format(' *** WARNING(111): File not found -- try again') goto 112 endif endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Output file? (Forced to 'yes' if there c - was an input file) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iinput.eq.1)then write(6,106) 106 format(/,1x,70('-'),/, * ' Write output to a file (y/n)? ',//, * ' -> ',$) read(5,'(a)')ofyn else ofyn = 'y' endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Open Output file if necessary ccccccccccccccccccccccccccccccccccccccccccccccccccc 1108 outfil = .false. if(ofyn.eq.'Y' .or. ofyn.eq.'y')then outfil = .true. write(6,107) 107 format(/,1x,70('-'),/, * ' Output file name?',//, * ' -> ',$) read(5,'(a)')ofnam open(lout,file=ofnam,status='new',form='formatted', * iostat = iss) c - File exists... if(iss.ne.0)then write(6,1107) 1107 format(/,' *** ERROR(1107): File Exists, try again',/) goto 1108 endif endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Begin interpolation procedure. Input file c - first, then keyboard ccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccc c - USING AN INPUT FILE ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iinput.eq.2)then ccccccccccccccccccccccccccccccccccccccccccccccccccc c - INPUT FILE, FREE FORMAT (FOR GEOID) TYPE 1 ccccccccccccccccccccccccccccccccccccccccccccccccccc if (iform.eq.1)then c - Read one line 113 format(a80) 116 read(linn,113,end=114)card kount = kount + 1 c - Find the lat/lon value c - Longitude always returns as 0->360...whether this c - is positive east or west is fixed a few lines down call ff1(card,xlat,xlon) c - If the lat/lon values came back as -999, set the c - geoid value to -999 and skip the interpolation if(xlat.eq.-999.d0 .or. xlon.eq.-999.d0)then val = -999.d0 goto 130 endif c - Force longitude to be positive East if(.not.poseast)xlon = 360.d0 - xlon c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Finally, write out the record to screen and possibly c - to an output file 130 call ff1out(iinput,outfil,lout,card,poseast,xlat,xlon,val) c - Go back and get another value goto 116 c - When finished, give a little report to the c - screen and end program. 114 ave = ave/keep rms = sqrt(rms/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) std = sqrt(fact*(rms**2 - ave**2)) else std = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * min,minlat,minlon,max,maxlat,maxlon, * ave,std,rms else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * min,minlat,360.d0 - minlon,max,maxlat,360.d0 - maxlon, * ave,std,rms endif stop 201 format(/,1x,70('-'),/, * ' FINAL REPORT: ',/, * ' Number of points input: ',i8,/, * ' Number of good points : ',i8,/, * ' Northernmost Latitude : ',f10.6,/, * ' Southernmost Latitude : ',f10.6,/, * ' Westernmost Longitude : ',f10.6,/, * ' Easternmost Longitude : ',f10.6,/, * ' Minimum Geoid Height : ',f8.3,/, * ' Lat/Lon of Minimum : ',f10.6,1x,f10.6,/, * ' Maximum Geoid Height : ',f8.3,/, * ' Lat/Lon of Maximum : ',f10.6,1x,f10.6,/, * ' Average Geoid Height : ',f8.3,/, * ' Standard Deviation : ',f8.3,/, * ' Root Mean Square : ',f8.3,//) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - INPUT FILE, FREE FORMAT (FOR GEOID) TYPE 2 ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(iform.eq.2)then c - Read one line 117 format(a80) 119 read(linn,117,end=118)card kount = kount + 1 c - Find the lat/lon value c - Longitude always returns as 0->360...whether this c - is positive east or west is fixed a few lines down call ff2(card,xlat,xlon) c - If the lat/lon values came back as -999, set the c - geoid value to -999 and skip the interpolation if(xlat.eq.-999.d0 .or. xlon.eq.-999.d0)then val = -999.d0 goto 330 endif c - Force longitude to be positive East if(.not.poseast)xlon = 360.d0 - xlon c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Finally, write out the record to screen and possibly c - to an output file 330 call ff2out(outfil,lout,card,poseast,xlat,xlon,val) c - Go back and get another value goto 119 c - When finished, give a little report to the c - screen and end program. 118 ave = ave/keep rms = sqrt(rms/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) std = sqrt(fact*(rms**2 - ave**2)) else std = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * min,minlat,minlon,max,maxlat,maxlon, * ave,std,rms else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * min,minlat,360.d0-minlon,max,maxlat,360.d0-maxlon, * ave,std,rms endif stop ccccccccccccccccccccccccccccccccccccccccccccccccccc c - INPUT FILE, BLUE BOOK FORMAT ccccccccccccccccccccccccccccccccccccccccccccccccccc else c - Spin through the blue book file, stopping at each c - *80* record, and while there, look forward c - to the next records to see if there is an *86* c - record. If so, overwrite the geoid field in c - it with an interpolated value. If not, c - create an *86* record. c - Read one line, and if it is NOT an *84* record, c - replicate it immediately into c - the output file (if we're using one). If the c - record is not an *80* record, come back and get c - a new record. 136 format(a80) 137 read(linn,136,end=138)card 153 if(card(8:9).eq.'84' .or. card(8:9).eq.'83')goto 137 if(outfil)then write(lout,136)card endif if(card(8:9).ne.'80')goto 137 c - arriving here, we've got an *80* record as 'card' c - Get the lat/lon value from the *80* record call bb80ll(card,xlat,xlon) c - If the xlat/xlon values are bad, don't do an interpolation. c - Just move on to the next record...any existing *86* c - record associated with this erroneous *80* record will c - remain unmodified. if(xlat.eq.-999.d0)goto 137 c - Arriving here, the *80* record has a good lat/lon value kount = kount + 1 c - The associated *86* record, if it exists, must c - come AFTER the *80* record we're currently c - looking at, and the ONLY thing allowed between the 80 c - and 86 record are 81, 82 and 85 records (83 and 84 c - records will be deleted). 140 read(linn,136,end=1139)card2 c - POSSIBILITY 1: We find an *83* or *84* record c - ACTION: Do not allow it to go into the output file, and c - go back for another "card2" value if(card2(8:9).eq.'83' .or. card2(8:9).eq.'84')then goto 140 endif c - POSSIBILITY 2: We find an *81*, *82* or *85* record c - ACTION: Replicate this record into the output file c - and go back for another "card2" value if(card2(8:9).eq.'81' .or. card2(8:9).eq.'82' .or. * card2(8:9).eq.'85')then if(outfil)then write(lout,136)card2 endif goto 140 endif c - POSSIBILITY 3: We find an *86* record c - ACTION: a) If it is the wrong *86* record, delete c it and make a NEW *86* record c and then go back and get another "card2" value c - b) Modify it if it's the right *86* record if(card2(8:9).eq.'86')then if(card2(11:14).eq.card(11:14))then c - Arriving here, we've found the *86* record associated with c - the *80* record c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Modify the *86* record and write it out c - Find the right geoid code c - GEOID99 if(imodel.eq.3)card2(43:43) = 'T' c - GSSS99 if(imodel.eq.4)card2(43:43) = 'U' c - GEOID03 if(imodel.eq.1)card2(43:43) = 'W' c - USGG2003 if(imodel.eq.2)card2(43:43) = 'X' ight = nint(val*1000) c write(card2(36:42),150)ight 150 format(i7) if(outfil)then write(lout,936)card2(1:35),ight,card2(43:80) endif c - And now go back and look for the next good *80* record in "card" goto 137 c - If this is NOT the associated *86* record, delete c - it and make a new *86* record else c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Make a NEW *86* record and write it out ight = nint(val*1000) card86 = b40//b40 card86(1:6) = card(1:6) card86(7:10) = '*86*' card86(11:14) = card(11:14) c - Find the right geoid code c - GEOID99 c2004apr16 if(imodel.eq.3) card2(43:43) = 'T' if(imodel.eq.3)card86(43:43) = 'T' c - GSSS99 c2004apr16 if(imodel.eq.4) card2(43:43) = 'U' if(imodel.eq.4)card86(43:43) = 'U' c - GEOID03 c2004apr16 if(imodel.eq.1) card2(43:43) = 'W' if(imodel.eq.1)card86(43:43) = 'W' c - USGG2003 c2004apr16 if(imodel.eq.2) card2(43:43) = 'X' if(imodel.eq.2)card86(43:43) = 'X' c write(card86(36:42),150)ight if(outfil)then write(lout,936)card86(1:35),ight,card86(43:80) endif c - And now go back and look for the next good *80* record in "card" goto 137 endif endif c - POSSIBILITY 4: We find something besides an 81,82, c - 83, 84, 85 or 86 record c - ACTION: Create a new *86* based on the previous *80* record, c - and then transfer the new 'card2' into 'card' and c - go back just after the 'read' statement for 'card' if(card2(8:9).ne.'81' .and. card2(8:9).ne.'82' .and. * card2(8:9).ne.'83' .and. card2(8:9).ne.'84' .and. * card2(8:9).ne.'85' .and. card2(8:9).ne.'86')then c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Make a NEW *86* record ight = nint(val*1000) card86 = b40//b40 card86(1:6) = card(1:6) card86(7:10) = '*86*' card86(11:14) = card(11:14) c write(card86(36:42),152)ight c 152 format(i7) c - Find the right geoid code c - GEOID99 c2004apr16 if(imodel.eq.3) card2(43:43) = 'T' if(imodel.eq.3)card86(43:43) = 'T' c - GSSS99 c2004apr16 if(imodel.eq.4) card2(43:43) = 'U' if(imodel.eq.4)card86(43:43) = 'U' c - GEOID03 c2004apr16 if(imodel.eq.1) card2(43:43) = 'W' if(imodel.eq.1)card86(43:43) = 'W' c - USGG2003 c2004apr16 if(imodel.eq.2) card2(43:43) = 'X' if(imodel.eq.2)card86(43:43) = 'X' c - ...and write out the *86* record if(outfil)then write(lout,936)card86(1:35),ight,card86(43:80) endif 936 format(a35,i7,a38) c - Now put the latest *80* record, which was in 'card2' c - into 'card' and go back to the top to search card = card2 goto 153 endif c - POSSIBILITY 5: We find the EOF c - ACTION: Create a new *86* based on the previous *80* record c - and then go to the 'end report' phase c - Find which geoid file to use, based on the lat/lon input 1139 call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Make a NEW *86* record ight = nint(val*1000) card86 = b40//b40 card86(1:6) = card(1:6) card86(7:10) = '*86*' card86(11:14) = card(11:14) c write(card86(36:42),152)ight c 152 format(i7) c - Find the right geoid code c - GEOID99 c2004apr16 if(imodel.eq.3) card2(43:43) = 'T' if(imodel.eq.3)card86(43:43) = 'T' c - GSSS99 c2004apr16 if(imodel.eq.4) card2(43:43) = 'U' if(imodel.eq.4)card86(43:43) = 'U' c - GEOID03 c2004apr16 if(imodel.eq.1) card2(43:43) = 'W' if(imodel.eq.1)card86(43:43) = 'W' c - USGG2003 c2004apr16 if(imodel.eq.2) card2(43:43) = 'X' if(imodel.eq.2)card86(43:43) = 'X' c - ...and write out the *86* record if(outfil)then write(lout,936)card86(1:35),ight,card86(43:80) endif c - Now go write out the report goto 138 c - This statement ought never to be reached, as all catches c - are made and returns coded with a catch (and, if you c - notice, upon compilation a "statement can't be reached" c - message appears) c goto 140 c - When finished, give a little report to the c - screen and end program. 138 ave = ave/keep rms = sqrt(rms/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) std = sqrt(fact*(rms**2 - ave**2)) else std = 0 endif write(6,201)kount,keep,xn,xs,xw,xe, * min,minlat,minlon,max,maxlat,maxlon, * ave,std,rms stop endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - USING KEYBOARD INPUT - Prompts ccccccccccccccccccccccccccccccccccccccccccccccccccc else c - Prompt user for station name 133 write(6,131) 131 format(/,1x,70('-'),/, * ' What is the name of this point?',//, * ' -> ',$) read(5,'(a)')sname kount = kount + 1 c - Prompt user for latitude 502 write(6,132) 132 format(/,1x,70('-'),/, * ' What is the latitude of this point?',//, * ' -> ',$) read(5,'(a)')xlatc call c2v(xlatc,xlat,1) if(xlat.eq.-999.d0)then write(6,501) 501 format(' *** ERROR(501): Bad Latitude...try again...') goto 502 endif c - Prompt user for longitude, using their desired convention 503 if(poseast)write(6,135) if(.not.poseast)write(6,1135) 135 format(/,1x,70('-'),/, * ' What is the longitude of this point?',/, * ' (Positive East, 0->360)',// * ' -> ',$) 1135 format(/,1x,70('-'),/, * ' What is the longitude of this point?',/, * ' (Positive West, 0->360)',// * ' -> ',$) c - Convert the card into variable "xlon" c - Longitude always returns as 0->360...whether this c - is positive east or west is fixed a few lines down read(5,'(a)')xlonc call c2v(xlonc,xlon,2) if(xlon.eq.-999.d0)then write(6,504) 504 format(' *** ERROR(504): Bad Longitude...try again...') goto 503 endif c - Force the longitude to be positive EAST if(.not.poseast)xlon = 360.d0 - xlon c - Find which geoid file to use, based on the lat/lon input call which1(xlat,xlon,nfiles,k,imodel) c - If the point isn't in any of our grid areas, set to -999 if(k.eq.-1)then val = -999 else c - Otherwise, do the interpolation call interp(xlat,xlon,lin,k,val) keep = keep + 1 ave = ave + val rms = rms + val*val if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon if(val.lt.min)then min = val minlat = xlat minlon = xlon endif if(val.gt.max)then max = val maxlat = xlat maxlon = xlon endif endif c - Finally, write out the record to screen and possibly c - to an output file. Use Free Format 1 output card = sname//b40 call ff1out(iinput,outfil,lout,card,poseast,xlat,xlon,val) c - Ask user if another point is to be used, and c - if so, go back and get it. write(6,134) 134 format(/,1x,70('-'),/, * ' Compute another?',// * ' -> ',$) read(5,'(a)')dumyn if(dumyn.eq.'y' .or. dumyn.eq.'Y')goto 133 c - When finished, give a little report to the c - screen and end program. 121 ave = ave/keep rms = sqrt(rms/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) std = sqrt(fact*(rms**2 - ave**2)) else std = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * min,minlat,minlon,max,maxlat,maxlon, * ave,std,rms else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * min,minlat,360.d0-minlon,max,maxlat,360.d0-maxlon, * ave,std,rms endif stop endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - End Of Main Program ccccccccccccccccccccccccccccccccccccccccccccccccccc end c c --------------------------------------------------------------- c subroutine files(imodel,dirnam,nfiles,fnam,lin,nff) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to create the file names for c - the gridded input data, and open all c - of those files that are available. ccccccccccccccccccccccccccccccccccccccccccccccccccc character*256 dirnam,b256 integer*4 lin(50) character*256 fnam(50) integer*4 ios(50) character*4 suf character*2 cval character*40 b40 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Use common statements to cut down on the c - amount of data sent through 'call' statements ccccccccccccccccccccccccccccccccccccccccccccccccccc common /grid3/ ios b40 = ' '// * ' ' b256=b40//b40//b40//b40//b40//b40//' ' if(dirnam.eq.b256)then ilen = 0 else ilen = lnblnk(dirnam) endif 109 format(i2.2) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Establish the suffix of files to be opened ccccccccccccccccccccccccccccccccccccccccccccccccccc suf = '.bin' cccccccccccccccccccccccccccccccccccccc c - The GEOID03 file names cccccccccccccccccccccccccccccccccccccc 2 if(imodel.eq.1)then nfiles = 14 c - First 8 files are CONUS do 301 i=1,8 write(cval,109)i fnam(i) = dirnam(1:ilen)//'g2003u'//cval//suf 301 continue cccccccccccccccccccccccccccccccccccccc c - The USGG2003 file names cccccccccccccccccccccccccccccccccccccc elseif(imodel.eq.2)then nfiles = 8 c - First 8 files are CONUS do 305 i=1,8 write(cval,109)i fnam(i) = dirnam(1:ilen)//'s2003u'//cval//suf 305 continue endif ccccccccccccccccccccccccccccccccccccc c - The GEOID99 file names cccccccccccccccccccccccccccccccccccccc if(imodel.eq.3)then nfiles = 14 c - First 8 files are CONUS do 101 i=1,8 write(cval,109)i fnam(i) = dirnam(1:ilen)//'g1999u'//cval//suf 101 continue c - Next 1 file is HAWAII do 102 i=9,9 write(cval,109)i-8 fnam(i) = dirnam(1:ilen)//'g1999h'//cval//suf 102 continue c - Next 1 file is PR/VI do 103 i=10,10 write(cval,109)i-9 fnam(i) = dirnam(1:ilen)//'g1999p'//cval//suf 103 continue c - Next 4 files are Alaska do 104 i=11,14 write(cval,109)i-10 fnam(i) = dirnam(1:ilen)//'g1999a'//cval//suf 104 continue cccccccccccccccccccccccccccccccccccccc c - The G99SSS file names cccccccccccccccccccccccccccccccccccccc elseif(imodel.eq.4)then nfiles = 8 c - First 8 files are CONUS do 105 i=1,8 write(cval,109)i fnam(i) = dirnam(1:ilen)//'s1999u'//cval//suf 105 continue cccccccccccccccccccccccccccccccccccccc c - Default Experimental Geoid: XXUSG cccccccccccccccccccccccccccccccccccccc elseif(imodel.eq.5)then nfiles = 8 c - First 8 files are CONUS do 777 i=1,8 write(cval,109)i fnam(i) = dirnam(1:ilen)//'xxusgu'//cval//suf 777 continue endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Now open all the files ccccccccccccccccccccccccccccccccccccccccccccccccccc do 201 i=1,nfiles open(lin(i),file=fnam(i),status='old',form='unformatted', * access='direct',recl=4,iostat=ios(i)) if(ios(i).eq.0)then write(6,203) write(6,204) fnam(i)(1:ilen+12) endif 201 continue 203 format(' *** Opening File: ',$) 204 format(a) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Check and see if at least ONE file was opened, c - and make a count of how many WERE opened. c - Abort if we find no geoid files. ccccccccccccccccccccccccccccccccccccccccccccccccccc nff = 0 do 202 i=1,nfiles if(ios(i).eq.0)nff = nff + 1 202 continue if(nff.gt.0)return write(6,205) 205 format(' *** ERROR(205): No files found -- aborting') stop end c c --------------------------------------------------------------- c subroutine expform ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to print an explanation of the various c - formats to the screen. ccccccccccccccccccccccccccccccccccccccccccccccccccc character*1 keyb ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Find out which format needs explaining ccccccccccccccccccccccccccccccccccccccccccccccccccc 2 write(6,1) 1 format(/,1x,70('-'),/, * ' Which format would you like explained? ',//, * ' 1 = Free Format (For Geoid) Type 1',/, * ' 2 = Free Format (For Geoid) Type 2',/, * ' 3 = NGS Blue Book Format',//, * ' 0 = End Program',/, * ' 99 = Return me to the main program',//, * ' -> ',$) read(5,*)iformx if(iformx.eq.0)stop if(iformx.eq.99)return ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Explain Free Format (For Geoid) Type 1 ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iformx.eq.1)then write(6,3) 3 format(/,1x,70('-'),/, * ' 1 = Free Format (For Geoid) Type 1 : ',/, * ' This is an ASCII format, where the first 40 characters',/, * ' of each record may contain the station name or be blank.',/, * ' The rest of the record (characters 41-80) must contain',/, * ' the latitude and longitude in one of three formats: ',/, * ' (a) decimal/integer degrees', * ' (two numbers)',/, * ' (b) integer degrees, decimal/integer minutes', * ' (four numbers)',/, * ' (c) integer degrees, integer minutes, decimal/integer', * ' seconds', * ' (six numbers)',/, * ' The latitude must be positive North. The longitude',/, * ' sign convention used is specified by the users during the',/, * ' running of this program.',// * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb write(6,4) 4 format(/,' EXAMPLES:',/ * ' The following records are examples of acceptable input:',//, * 1x,40x,'<------------ Columns 41-80 ----------->' ,/, * ' BM1027',34x,'50 285',/, * ' XX8063',34x,' 45.3 282',/, * 1x,16x,'AIRPORT SC1',13x,' 47 18.55 252 14',/, * 1x,40x,'30 17 270 59',/, * 1x,10x,'CORNER MARKER',17x,'48 14 17 239 52 18.753',/, * ' DATA PT #4077',27x,' 42 42 42.4242 255 55 55.55555',// * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb goto 2 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Explain Free Format (For Geoid) Type 2 ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iformx.eq.2)then write(6,5) 5 format(/,1x,70('-'),/, * ' 2 = Free Format (For Geoid) Type 2 : ',/, * ' This is an ASCII format, where the last 40 characters',/, * ' of each record may contain the station name or be blank.',/, * ' The first 32 characters of the record (characters 1-32)',/, * ' must contain the latitude and longitude in one of three',/, * ' formats:',/, * ' (a) decimal/integer degrees', * ' (two numbers)',/, * ' (b) integer degrees, decimal/integer minutes', * ' (four numbers)',/, * ' (c) integer degrees, integer minutes, decimal/integer', * ' seconds', * ' (six numbers)',/, * ' The latitude must be positive North. The longitude',/, * ' sign convention used is specified by the users during the',/, * ' running of this program.',//, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb write(6,6) 6 format(/,' EXAMPLES:',/ * ' The following records are examples of acceptable input:',//, * ' <-------- Columns 1-32 -------->',/, * ' 50 285',26x,'BM1027',/, * ' 45.3 282',15x,'XXX8063',/, * ' 47 18.55 252 14',31x,'AIRPORT SC1',/, * ' 30 17 270 59',/, * ' 48 14 17 239 52 18.753',20x,'CORNER MARKER',/, * ' 42 42 42.4242 255 55 55.55555',1x,'DATA PT #4077',//, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb goto 2 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Explain NGS Blue Book Format ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iformx.eq.3)then write(6,7) 7 format(/,1x,70('-'),/, * ' 3 = NGS Blue Book Format : ',/) WRITE (6,510) 510 FORMAT (' NGS Blue Book format - *80*', + ' (Control Point) and *86* (Combo', /, + ' Height) Records. INTG uses information from the', + ' *80* and *86* records in Blue', /, + ' Book files. The Station Serial Number (SSN), the', + ' latitude and the longitude', /, + ' are read from the *80* records, the rest of the', + ' record is ignored. Only the', /, + ' *86* record is modified by INTG - the *80*', + ' records and all other records are') WRITE (6, 511) 511 FORMAT (' passed through without change to the output file.', + ' Any *84* records are ', /, + ' deleted. Information in columns', + ' 36-43 in the *86* record is', /, + ' created by INTG. ',/) WRITE (6, 520) 520 FORMAT (' For more information on this format,', + ' please refer to:', /, + ' ''Input Formats and Specifications of the', + ' National Geodetic Survey Data Base''', /, + ' ''Volume 1. Horizontal Control Data''.', /, + ' Published by the Federal Geodetic Control Committee', + ' in September 1994 and', /, + ' available from: the National Geodetic Survey, NOAA,', + ' Silver Spring, MD 20910.', //, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb WRITE (6, 530) 530 FORMAT (' The following input example is *80* and *86* records', + ' from a Blue Book', /, + ' file:', //, + ' 004560*80*0096KNOXVILLE CT HSE ', + '411906578 N0930548534 W 277 MIA33', /, + ' 004565*86*0096 any geoid values in here ', + 'will be overwritten if ssn matches ',//, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb goto 2 endif end c c --------------------------------------------------------------- c subroutine intro ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to print out introductory screens c - and disclaimers ccccccccccccccccccccccccccccccccccccccccccccccccccc character*1 keyb ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Introduction ccccccccccccccccccccccccccccccccccccccccccccccccccc write(6,1) 1 format(////////////, * 10x,' Welcome to the ',/, * 10x,' National Geodetic Survey''s ',/, * 10x,' INTG PROGRAM',/, * 10x,' (Interpolate from Geoid files). ',/) write(6,2) 2 format( * 10x,' VERSION 1.3',/, * 10x,' August 29, 2003',/, * 10x,' Daniel R. Roman, Ph.D.',//, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Disclaimer ccccccccccccccccccccccccccccccccccccccccccccccccccc WRITE (6,932) 932 format(/,1x,70('-'),/, * /, 32X, 'DISCLAIMER' ,//, + ' This program and supporting information is furnished by', + ' the government of', /, + ' the United States of America, and is accepted/used by the', + ' recipient with', /, + ' the understanding that the U. S. government makes no', + ' warranties, express or', /, + ' implied, concerning the accuracy, completeness, reliability,', + ' or suitability', /, + ' of this program, of its constituent parts, or of any', + ' supporting data.', //, + ' The government of the United States of America shall be', + ' under no liability', /, + ' whatsoever resulting from any use of this program.', + ' This program should', /, + ' not be relied upon as the sole basis for solving a problem', + ' whose incorrect', /, + ' solution could result in injury to person or property.') WRITE (6,933) 933 FORMAT ( /, + ' This program is the property of the government of the', + ' United States of', /, + ' America. Therefore, the recipient further agrees not to', + ' assert proprietary', /, + ' rights therein and not to represent this program to anyone as', + ' being other', /, + ' than a government program.', //, * ' (Hit RETURN to continue)') read(5,'(a)')keyb WRITE (6,934) 934 FORMAT ( /, + ' The flag for bad data is ''-999'' which will mostly appear if', + ' the input ', /, + ' latitude/longitude of a point falls outside of the borders', + ' of available files.', /, + ' Note the longitude can be entered in either positive EAST', + ' or ', /, + ' as positive WEST.',//, * ' Allowable formats for latitude and longitude, both keyboard', * ' and by file are:',/, * ' 1) Integer degrees 2) Decimal degrees',/, * ' 3) Int Deg & Integer Minutes 4) Int Deg & Decimal', * ' Minutes',/, * ' 5) Int Deg, Int Min, & Integer Sec 6) Int Deg, Int Min &', * ' Decimal Seconds',/, * ' Space delimited values only. No commas.',/, * ' Negative values are allowed, but only with the negative ', * ' sign preceding',/, * ' and touching degrees (i.e. -100 42 37.3)',//, * ' (Hit RETURN to continue)') read(5,'(a)')keyb return end c c --------------------------------------------------------------- c subroutine ff1(card,xlat,xlon) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to take one "Free Format, Type 1" c - record (card) and extract the latitude and c - longitude (xlat,xlon) in decimal degrees ccccccccccccccccccccccccccccccccccccccccccccccccccc character*80 card character*40 ecard real*8 xlat,xlon integer*4 bkt,ekt integer*4 b(50),e(50) real*8 deglat,minlat,seclat real*8 deglon,minlon,seclon integer*4 ival,iflag real*8 xval bkt = 0 ekt = 0 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Search for the beginnings c - and endings of numbers, assuming only c - that we'll find the numbers '0' through '9' as c - well as spaces, ' ', and decimals '.', and c - perhaps a leading negative sign. c - c - Comma delimeted data will not work ccccccccccccccccccccccccccccccccccccccccccccccccccc ecard = card(41:80) ilen = lnblnk(ecard) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - The position of the character which begins c - the first 'number' is b(1), and the position c - of the character which ends the first 'number' c - is e(1).....etc for all numbers. (Positions c - are counted, in Free Format Type 1, as counting c - 1 to 40(at most), starting in the 41st character c - of the record.) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(ecard(1:1).ne.' ')then bkt = bkt + 1 b(bkt) = 1 endif do 1 i=2,ilen if(ecard(i:i).eq.' ')then if(ecard(i-1:i-1).ne.' ')then ekt = ekt + 1 e(ekt) = i-1 endif else if(ecard(i-1:i-1).eq.' ')then bkt = bkt + 1 b(bkt) = i endif endif 1 continue c - Count the last space as an end, since we've c - already defined its length as the last c - non-space character ekt = ekt + 1 e(ekt) = i-1 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Obviously the number of beginnings and c - endings must be the same. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .ne. ekt)stop 80808 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer/Decimal Degrees (2 numbers) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .eq. 2)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,deglon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval xlat = deglat xlon = deglon ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees & Integer/Decimal c - minutes (4 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 4)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,minlat,deglon,minlon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)minlat = ival if(iflag.eq.1)minlat = xval if(deglat.lt.0)minlat = -minlat call val(ecard(b(3):e(3)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval call val(ecard(b(4):e(4)),iflag,ival,xval) if(iflag.eq.0)minlon = ival if(iflag.eq.1)minlon = xval if(deglon.lt.0)minlon = -minlon xlat = deglat+(minlat/60.d0) xlon = deglon+(minlon/60.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees, Integer Minutes & c - Integer/Decimal Seconds (6 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 6)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,minlat,seclat, c * deglon,minlon,seclon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)minlat = ival if(iflag.eq.1)minlat = xval if(deglat.lt.0)minlat = -minlat call val(ecard(b(3):e(3)),iflag,ival,xval) if(iflag.eq.0)seclat = ival if(iflag.eq.1)seclat = xval if(deglat.lt.0)seclat = -seclat call val(ecard(b(4):e(4)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval call val(ecard(b(5):e(5)),iflag,ival,xval) if(iflag.eq.0)minlon = ival if(iflag.eq.1)minlon = xval if(deglon.lt.0)minlon = -minlon call val(ecard(b(6):e(6)),iflag,ival,xval) if(iflag.eq.0)seclon = ival if(iflag.eq.1)seclon = xval if(deglon.lt.0)seclon = -seclon xlat = deglat+(minlat/60.d0)+(seclat/3600.d0) xlon = deglon+(minlon/60.d0)+(seclon/3600.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Otherwise there's an error and the lat/lon c - are set to -999 ccccccccccccccccccccccccccccccccccccccccccccccccccc else xlat = -999.d0 xlon = -999.d0 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If, somehow, the lat value isn't in the c - right range, turn both lat/lon into -999s. c - If, somehow, the lon value isn't in the c - right range, try some simple arithmetic c - and if that fails, turn both lat/lon into c - -999s ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlon.lt.000)xlon = xlon + 360.d0 if(xlon.gt.360)xlon = xlon - 360.d0 if(xlat.lt.-90 .or. xlat.gt.+90 .or. * xlon.lt.000 .or. xlon.gt.360)then xlat = -999.d0 xlon = -999.d0 endif return end c c --------------------------------------------------------------- c subroutine ff2(card,xlat,xlon) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to take one "Free Format, Type 2" c - record (card) and extract the latitude and c - longitude (xlat,xlon) in decimal degrees ccccccccccccccccccccccccccccccccccccccccccccccccccc character*80 card character*32 ecard real*8 xlat,xlon integer*4 bkt,ekt integer*4 b(50),e(50) real*8 deglat,minlat,seclat real*8 deglon,minlon,seclon integer*4 ival,iflag real*8 xval bkt = 0 ekt = 0 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Search for the beginnings c - and endings of numbers, assuming only c - that we'll find the numbers '0' through '9' as c - well as spaces, ' ', and decim, and c - perhaps a leading negative sign. c - c - Comma delimeted data will not work ccccccccccccccccccccccccccccccccccccccccccccccccccc ecard = card(1:32) ilen = lnblnk(ecard) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - The position of the character which begins c - the first 'number' is b(1), and the position c - of the character which ends the first 'number' c - is e(1).....etc for all numbers. (Positions c - are counted, in Free Format Type 2, as counting c - 1 to 32(at most), starting in the 1st character c - of the record.) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(ecard(1:1).ne.' ')then bkt = bkt + 1 b(bkt) = 1 endif do 1 i=2,ilen if(ecard(i:i).eq.' ')then if(ecard(i-1:i-1).ne.' ')then ekt = ekt + 1 e(ekt) = i-1 endif else if(ecard(i-1:i-1).eq.' ')then bkt = bkt + 1 b(bkt) = i endif endif 1 continue c - Count the last space as an end, since we've c - already defined its length as the last c - non-space character ekt = ekt + 1 e(ekt) = i-1 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Obviously the number of beginnings and c - endings must be the same. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .ne. ekt)stop 80808 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer/Decimal Degrees (2 numbers) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .eq. 2)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,deglon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval xlat = deglat xlon = deglon ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees & Integer/Decimal c - minutes (4 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 4)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,minlat,deglon,minlon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)minlat = ival if(iflag.eq.1)minlat = xval if(deglat.lt.0)minlat = -minlat call val(ecard(b(3):e(3)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval call val(ecard(b(4):e(4)),iflag,ival,xval) if(iflag.eq.0)minlon = ival if(iflag.eq.1)minlon = xval if(deglon.lt.0)minlon = -minlon xlat = deglat+(minlat/60.d0) xlon = deglon+(minlon/60.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees, Integer Minutes & c - Integer/Decimal Seconds (6 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 6)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deglat,minlat,seclat, c * deglon,minlon,seclon call val(ecard(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deglat = ival if(iflag.eq.1)deglat = xval call val(ecard(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)minlat = ival if(iflag.eq.1)minlat = xval if(deglat.lt.0)minlat = -minlat call val(ecard(b(3):e(3)),iflag,ival,xval) if(iflag.eq.0)seclat = ival if(iflag.eq.1)seclat = xval if(deglat.lt.0)seclat = -seclat call val(ecard(b(4):e(4)),iflag,ival,xval) if(iflag.eq.0)deglon = ival if(iflag.eq.1)deglon = xval call val(ecard(b(5):e(5)),iflag,ival,xval) if(iflag.eq.0)minlon = ival if(iflag.eq.1)minlon = xval if(deglon.lt.0)minlon = -minlon call val(ecard(b(6):e(6)),iflag,ival,xval) if(iflag.eq.0)seclon = ival if(iflag.eq.1)seclon = xval if(deglon.lt.0)seclon = -seclon xlat = deglat+(minlat/60.d0)+(seclat/3600.d0) xlon = deglon+(minlon/60.d0)+(seclon/3600.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Otherwise there's an error and the lat/lon c - must be re-input ccccccccccccccccccccccccccccccccccccccccccccccccccc else xlat = -999.d0 xlon = -999.d0 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If, somehow, the lat value isn't in the c - right range, turn both lat/lon into -999s. c - If, somehow, the lon value isn't in the c - right range, try some simple arithmetic c - and if that fails, turn both lat/lon into c - -999s ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlon.lt.000)xlon = xlon + 360.d0 if(xlon.gt.360)xlon = xlon - 360.d0 if(xlat.lt.-90 .or. xlat.gt.+90 .or. * xlon.lt.000 .or. xlon.gt.360)then xlat = -999.d0 xlon = -999.d0 endif return end c c --------------------------------------------------------------- c subroutine which1(xlat,xlon,nfiles,k,imodel) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to decide which of the open c - geoid files will be used to interpolate c - to a point at xlat/xlon. The returned c - value of "k" means the "kth" file will c - be used. c - c - For the GEOID99 models: c Alaska and CONUS overlap, and this code c - forces a "CONUS wins" scenario when interpolating c - the geoid at points in the overlap region. ccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 xlat,xlon integer*4 nfiles,k,imodel real*8 glamn(50),glamx(50),glomn(50),glomx(50) real*8 dla(50),dlo(50) integer*4 nla(50),nlo(50),ikind(50) integer*4 ios(50) integer*4 rank(50) logical ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Use common statements to cut down on the c - amount of data sent through 'call' statements ccccccccccccccccccccccccccccccccccccccccccccccccccc common /grid1/ glamn,glamx,glomn,glomx common /grid2/ dla,dlo,nla,nlo,ikind common /grid3/ ios common /edges/ ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Spin through all *open* files, and c - *RANK* them...the file with the c - highest RANK decides the value of "k". c - Here's how the ranking goes: c - 0 = Point does not lie in this file area c - 1 = Point lies in this file, but at a corner c - 2 = Point lies in this file, but at an edge c - 3 = Point lies in this file, away from corners/edges c - If a rank=3 file is found, k is immediately set c - and we return. If no 3 files are found, the c - file with the highest rank (2 or 1) sets k. c - If all files have rank=0, k is set to -1 ccccccccccccccccccccccccccccccccccccccccccccccccccc do 1 i=1,nfiles if(ios(i).ne.0)goto 1 rank(i) = 0 if(xlat.le.glamx(i) .and. xlat.ge.glamn(i) .and. * xlon.le.glomx(i) .and. xlon.ge.glomn(i) )then c - At this point, we're Inside a grid c - Now determine which one of the 9 possible c - places this point resides -- c - NW corner, North Edge, NE corner c - West Edge, Center , East edge c - SW corner, South Edge, SE corner ne = .false. se = .false. we = .false. ee = .false. c - Near North edge? if(glamx(i) - xlat .le. dla(i)/2)then ne=.true. c - Near South edge? elseif(xlat - glamn(i) .le. dla(i)/2)then se=.true. endif c - Near East edge? if(glomx(i) - xlon .le. dlo(i)/2)then ee=.true. c - Near West edge? elseif(xlon - glomn(i) .le. dlo(i)/2)then we=.true. endif c - If we're not near an edge, set "k" to "i" and return c - (The only exception to this is if we're using GEOID99 c - and have both Alaska and USA grids available and c - we're in the overlap region) c - This is the Alaska/CONUS overlap exception code if(imodel.eq.1)then if(i.ne.1 .and. ios(1).eq.0 .and. xlat.le.58 .and. * xlat.ge.49 .and. xlon.ge.230 .and. xlon.le.234)then k = 1 return endif endif if(.not.ne .and. .not.se .and. .not.we .and. .not.ee)then k = i return endif c - Set the rank of this file, based on edge-logic if(ne .and. .not.we .and. .not.ee)rank(i) = 2 if(se .and. .not.we .and. .not.ee)rank(i) = 2 if(we .and. .not.ne .and. .not.se)rank(i) = 2 if(ee .and. .not.ne .and. .not.se)rank(i) = 2 if(ne.and.we)rank(i) = 1 if(se.and.we)rank(i) = 1 if(se.and.ee)rank(i) = 1 if(ne.and.ee)rank(i) = 1 endif 1 continue c - If we reach this point, all possible files have c - been searched, and there's no open file which c - had a rank of 3. So now, see if we have any rank 2 c - or rank 1 files to use. do 101 i=1,nfiles if(ios(i).ne.0)goto 101 if(rank(i).eq.2)then k = i return endif 101 continue do 102 i=1,nfiles if(ios(i).ne.0)goto 102 if(rank(i).eq.1)then k = i return endif 102 continue c - If we come here, no files are acceptable for the c - given lat/lon value k=-1 return end c c --------------------------------------------------------------- c subroutine ff1out(iinput,outfil,lout,card,poseast,xlat,xlon,val) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to write out a "Free Format, (For c - Geoid) Type 1" record ccccccccccccccccccccccccccccccccccccccccccccccccccc logical outfil logical poseast integer*4 lout,iinput character*80 card character*40 fcard real*8 xlat,xlon real*4 val integer*4 deglat,minlat,deglon,minlon real*4 seclat,seclon fcard = card(1:40) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Compute Degrees/Minutes/Seconds from Lat/Lon c - (which is the default output for FF1) c - Make sure to convert to Positive WEST if that c - is what the user has chosen to work in. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlat.ne.-999 .and. xlon .ne. -999)then if(.not.poseast)xlon = 360.d0 - xlon deglat = int(xlat + 1.d-12) deglon = int(xlon + 1.d-12) minlat = int((xlat-dble(deglat))*60.d0 + 1.d-9) minlon = int((xlon-dble(deglon))*60.d0 + 1.d-9) seclat = (((xlat-deglat)*60.d0)-minlat)*60.d0 seclon = (((xlon-deglon)*60.d0)-minlon)*60.d0 else deglat = 99 minlat = 99 seclat = 99.9999999 deglon = 99 minlon = 99 seclon = 99.9999999 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If we are outputting to an output file, c - do so here ccccccccccccccccccccccccccccccccccccccccccccccccccc if(outfil)then write(lout,99)fcard,deglat,minlat,seclat, * deglon,minlon,seclon,val endif if(iinput.eq.1)write(6,98)fcard,deglat,minlat,seclat, * deglon,minlon,seclon,val 98 format(/,' Your Result:',/, * 1x,40x,' latitude ',1x, * ' longitude ',1x,' N ',/, * ' Station Name',26x,'ddd mm ss.sssss',1x, * 'ddd mm ss.sssss'1x,' meters',/, * 1x,a40,i3,1x,i2,1x,f8.5,1x, * i3,1x,i2,1x,f8.5,1x,f8.3,/) 99 format(a40,i3,1x,i2,1x,f8.5,1x, * i3,1x,i2,1x,f8.5,1x,f8.3) return end c c --------------------------------------------------------------- c subroutine ff2out(outfil,lout,card,poseast,xlat,xlon,val) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to write out a "Free Format, (For c - Geoid) Type 2" record ccccccccccccccccccccccccccccccccccccccccccccccccccc logical outfil logical poseast integer*4 lout character*80 card character*40 ecard real*8 xlat,xlon real*4 val integer*4 deglat,minlat,deglon,minlon real*4 seclat,seclon ecard = card(41:80) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Compute Degrees/Minutes/Seconds from Lat/Lon c - (which is the default output for FF2) c - Make sure to convert to Positive WEST if that c - is what the user has chosen to work in. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlat.ne.-999 .and. xlon .ne. -999)then if(.not.poseast)xlon = 360.d0 - xlon deglat = int(xlat + 1.d-12) deglon = int(xlon + 1.d-12) minlat = int((xlat-dble(deglat))*60.d0 + 1.d-9) minlon = int((xlon-dble(deglon))*60.d0 + 1.d-9) seclat = (((xlat-deglat)*60.d0)-minlat)*60.d0 seclon = (((xlon-deglon)*60.d0)-minlon)*60.d0 else deglat = 99 minlat = 99 seclat = 99.9999999 deglon = 99 minlon = 99 seclon = 99.9999999 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If we are outputting to an output file, c - do so here ccccccccccccccccccccccccccccccccccccccccccccccccccc if(outfil)then write(lout,99)deglat,minlat,seclat, * deglon,minlon,seclon,val,ecard endif 99 format(i3,1x,i2,1x,f8.5,1x, * i3,1x,i2,1x,f8.5,1x,f8.3,a40) return end c c --------------------------------------------------------------- c subroutine c2v(card,value,ilatlon) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine for converting a 40 character c - string into a value. The string must c - consist of 1,2 or 3 numbers, meaning c - Degrees, Degrees&Minutes or Deg, Min & Sec c - c - This has nothing to do with Positive c - Longitude conventions. It is purely a c - character-to-value conversion. Longitude c - conventions are handled outside this c - subroutine. c - c - The allowable characters in "card" are c - spaces, numbers 1 through 9, decimal c - points, and one leading negative sign. c - c - Returned values of "value" will be forced into: c - ilatlon = 1 => -90 < value < +90 c - ilatlon = 2 => 0 < value < +360 ccccccccccccccccccccccccccccccccccccccccccccccccccc character*40 card integer*4 ilatlon integer*4 bkt,ekt integer*4 b(50),e(50) real*8 deg,min,sec real*8 value integer*4 ival,iflag real*8 xval bkt = 0 ekt = 0 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Search for the beginnings c - and endings of numbers, assuming only c - that we'll find the numbers '0' through '9' as c - well as spaces, ' ', and decimals '.', and c - possibly a leading negative sign '-'. c - c - Comma delimeted data will not work c - A space between the negative sign and the c - first number is not allowable. ccccccccccccccccccccccccccccccccccccccccccccccccccc ilen = lnblnk(card) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - The position of the character which begins c - the first 'number' is b(1), and the position c - of the character which ends the first 'number' c - is e(1).....etc for all numbers. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(card(1:1).ne.' ')then bkt = bkt + 1 b(bkt) = 1 endif do 1 i=2,ilen if(card(i:i).eq.' ')then if(card(i-1:i-1).ne.' ')then ekt = ekt + 1 e(ekt) = i-1 endif else if(card(i-1:i-1).eq.' ')then bkt = bkt + 1 b(bkt) = i endif endif 1 continue c - Count the last space as an end, since we've c - already defined its length as the last c - non-space character ekt = ekt + 1 e(ekt) = i-1 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Obviously the number of beginnings and c - endings must be the same. ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .ne. ekt)stop 90909 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer/Decimal Degrees (1 number) ccccccccccccccccccccccccccccccccccccccccccccccccccc if(bkt .eq. 1)then c - The following line doesn't work with older FORTRAN compilers c read(card,*)deg call val(card(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deg = ival if(iflag.eq.1)deg = xval value = deg ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees & Integer/Decimal c - minutes (2 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 2)then c - The following line doesn't work with older FORTRAN compilers c read(card,*)deg,min call val(card(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deg = ival if(iflag.eq.1)deg = xval call val(card(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)min = ival if(iflag.eq.1)min = xval if(deg.lt.0)min = -min value = deg+(min/60.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Integer Degrees, Integer Minutes & c - Integer/Decimal Seconds (3 numbers) c - Take care with pos/neg signs ccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(bkt .eq. 3)then c - The following line doesn't work with older FORTRAN compilers c read(ecard,*)deg,min,sec call val(card(b(1):e(1)),iflag,ival,xval) if(iflag.eq.0)deg = ival if(iflag.eq.1)deg = xval call val(card(b(2):e(2)),iflag,ival,xval) if(iflag.eq.0)min = ival if(iflag.eq.1)min = xval if(deg.lt.0)min = -min call val(card(b(3):e(3)),iflag,ival,xval) if(iflag.eq.0)sec = ival if(iflag.eq.1)sec = xval if(deg.lt.0)sec = -sec value = deg+(min/60.d0)+(sec/3600.d0) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Otherwise there's an error and the value c - must be re-input ccccccccccccccccccccccccccccccccccccccccccccccccccc else value = -999.d0 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - If, somehow, the latitude isn't in the c - right range, turn lat/lon into -999s. c - If somehow, the longitude isn't in the c - right range, try a quick arithmetic trick, c - and then if it's still bad, turn the c - lat/lon values into -999s ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Latitude if(ilatlon.eq.1)then if(value.lt.-90 .or. value.gt.+90)value = -999.d0 c - Longitude else if(value.lt.000)value = value + 360.d0 if(value.gt.360)value = value - 360.d0 if(value.lt.000 .and. value.gt.360)value = -999.d0 endif return end c c --------------------------------------------------------------- c subroutine bb80ll(card,xlat,xlon) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine for extracting the latitude and c - longitude from a *80* blue book record c - Force: -90 <= xlat <= +90 c - 0 <= xlon <= +360 ccccccccccccccccccccccccccccccccccccccccccccccccccc character*80 card integer*4 deglat,minlat,seclat integer*4 deglon,minlon,seclon character*1 lathem,lonhem real*8 xlat,xlon,latsc,lonsc integer*4 ival,iflag real*8 xval c - The following line doesn't work with older FORTRAN compilers c read(card,1)deglat,minlat,seclat,lathem, c * deglon,minlon,seclon,lonhem 1 format(44x,i2,i2,i7,a1,i3,i2,i7,a1) call val(card(45:46),iflag,ival,xval) deglat = ival call val(card(47:48),iflag,ival,xval) minlat = ival call val(card(49:55),iflag,ival,xval) seclat = ival lathem = card(56:56) call val(card(57:59),iflag,ival,xval) deglon = ival call val(card(60:61),iflag,ival,xval) minlon = ival call val(card(62:68),iflag,ival,xval) seclon = ival lonhem = card(69:69) if(lathem.eq.'N')then latsc = +1 elseif(lathem.eq.'S')then latsc = -1 else xlat = -999.d0 xlon = -999.d0 return endif if(lonhem.eq.'E')then lonsc = +1 elseif(lonhem.eq.'W')then lonsc = -1 else xlat = -999.d0 xlon = -999.d0 return endif xlat = latsc*(deglat + (minlat/60.d0) + (seclat/3600.d0/1d5)) xlon = lonsc*(deglon + (minlon/60.d0) + (seclon/3600.d0/1d5)) if(xlon.lt.0)xlon = xlon + 360.d0 if(xlat.lt.-90 .or. xlat.gt.+90 .or. * xlon.lt.000 .or. xlon.gt.360) then xlat = -999.d0 xlon = -999.d0 endif return end c c c subroutine interp(xlat,xlon,lin,k,val) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine for interpolating a value(val) c - off of the kth grid. As the gridded data c - are all direct-access files, only the c - nearest points (1 to 9, depending on c - the points location relative to corners c - and edges) are read into RAM for each c - point's interpolation. c - c - The size/spacing/etc of the data grid c - are defined by the common statement, and c - the variable "k" c - Caveats: c It is assumed that xlat/xlon fall in the c region of the kth grid ccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 xlat,xlon integer*4 lin(50) integer*4 k real*4 val real*8 glamn(50),glamx(50),glomn(50),glomx(50) real*8 dla(50),dlo(50) integer*4 nla(50),nlo(50) integer*4 ikind(50) real*4 f1,f2,f3,f4,f5,f6,f7,f8,f9 logical ne,se,we,ee c variables and functions for handling endian problem logical reverse_bytes(50) common /endian/reverse_bytes real*4 frev4 common /grid1/ glamn,glamx,glomn,glomx common /grid2/ dla,dlo,nla,nlo,ikind common /edges/ ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - A little check for safety's sake ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlat.gt.glamx(k) .or. xlat.lt.glamn(k) .or. * xlon.lt.glomn(k) .or. xlon.gt.glomx(k))then stop 10000 endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Find the row/col of the nearest point to c - this lat/lon point ccccccccccccccccccccccccccccccccccccccccccccccccccc irown = nint((xlat-glamn(k)) / dla(k)) + 1 icoln = nint((xlon-glomn(k)) / dlo(k)) + 1 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Find the latitude and longitude of the nearest c - grid point ccccccccccccccccccccccccccccccccccccccccccccccccccc xlatn = glamn(k) + (irown-1)*dla(k) xlonn = glomn(k) + (icoln-1)*dlo(k) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - First things first -- if we're sitting right c - on a grid node, just assign the value and c - return ccccccccccccccccccccccccccccccccccccccccccccccccccc if(xlat.eq.xlatn .and. xlon.eq.xlonn)then irec = 11 + (irown-1)*nlo(k) + icoln read(lin(k),rec=irec)f1 c read big endian values and convert to little endian if(reverse_bytes(k)) f1=frev4(f1) val = f1 return endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Check and see if we're NEAR or ON an edge... c - c - Although we'll stick with biquadratic c - interpolation, we've got to pick 9 points c - that are inside the grid, and therefore are c - not the exact 9 that would have been picked c - if the grid edge weren't hindering us. c - But there's no choice, since the program c - logic forces an edge check before coming c - to this subroutine....i.e., if we're here, c - and near an edge, there's no data outside c - of the edge we're near. c - c - If we ARE near or on an edge, then c - the values "irown" and "icoln" will now c - be changed to no longer mean "nearest" c - grid node, but will now mean "central node" c - of the 9x9 points we MUST pick for doing c - near-edge computations ccccccccccccccccccccccccccccccccccccccccccccccccccc if(irown.eq.1) irown = 2 if(irown.eq.nla(k))irown = nla(k) - 1 if(icoln.eq.1) icoln = 2 if(icoln.eq.nlo(k))icoln = nlo(k) - 1 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - At this point, whether we're on an edge or c - not, the the "irown/icoln" values reflect c - the center node of the 9x9 points we have to c - use for the biquadratic interpolation ccccccccccccccccccccccccccccccccccccccccccccccccccc irec = 11 + (irown-1-1)*nlo(k) + icoln - 1 read(lin(k),rec=irec)f1 irec = 11 + (irown-1-1)*nlo(k) + icoln read(lin(k),rec=irec)f2 irec = 11 + (irown-1-1)*nlo(k) + icoln + 1 read(lin(k),rec=irec)f3 irec = 11 + (irown -1)*nlo(k) + icoln - 1 read(lin(k),rec=irec)f4 irec = 11 + (irown -1)*nlo(k) + icoln read(lin(k),rec=irec)f5 irec = 11 + (irown -1)*nlo(k) + icoln + 1 read(lin(k),rec=irec)f6 irec = 11 + (irown+1-1)*nlo(k) + icoln - 1 read(lin(k),rec=irec)f7 irec = 11 + (irown+1-1)*nlo(k) + icoln read(lin(k),rec=irec)f8 irec = 11 + (irown+1-1)*nlo(k) + icoln + 1 read(lin(k),rec=irec)f9 c convert big endian values to little endian if(reverse_bytes(k)) then f1=frev4(f1) f2=frev4(f2) f3=frev4(f3) f4=frev4(f4) f5=frev4(f5) f6=frev4(f6) f7=frev4(f7) f8=frev4(f8) f9=frev4(f9) endif c f1 = data(irown-1 , icoln-1) c f2 = data(irown-1 , icoln ) c f3 = data(irown-1 , icoln+1) c f4 = data(irown , icoln-1) c f5 = data(irown , icoln ) c f6 = data(irown , icoln+1) c f7 = data(irown+1 , icoln-1) c f8 = data(irown+1 , icoln ) c f9 = data(irown+1 , icoln+1) xx = ( xlon - (glomn(k)+(icoln-2)*dlo(k)) ) / dlo(k) yy = ( xlat - (glamn(k)+(irown-2)*dla(k)) ) / dla(k) fx1 = qfit(xx,f1,f2,f3) fx2 = qfit(xx,f4,f5,f6) fx3 = qfit(xx,f7,f8,f9) val = qfit(yy,fx1,fx2,fx3) return end c c --------------------------------------------------------------- c real function qfit(x,f0,f1,f2) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Parabola fit through three points (x=0, x=1, x=2) c - with values f0=f(0), f1=f(1), and f2=f(2) c - and returning the value c - qfit = f(x), where 0<=x<=2 ccccccccccccccccccccccccccccccccccccccccccccccccccc df0 = f1 - f0 df1 = f2 - f1 d2f0 = df1 - df0 qfit = f0 + x*df0 + 0.5*x*(x-1)*d2f0 return end c c --------------------------------------------------------------- c integer function lnblnk(card) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Function to return the position of the last c - non - blank character of a string. c - This function is found in most FORTRAN c - language compilers. It is included here c - because some compilers do not yet recognize c - it. ccccccccccccccccccccccccccccccccccccccccccccccccccc character *(*) card il = len(card) do 1 i=1,il if(card(i:i).ne.' ')ix = i 1 continue lnblnk = ix return end c c --------------------------------------------------------------- c subroutine val(card,it,iv,xv) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to return the value of a character c - string as either an integer (iv) or a c - real (xv), with "it" telling which value c - to use (it = 0 means int, it=1 means real) c - This subroutine is included to alleviate c - the trouble that old compilers have in c - doing simple character manipulations c - For now, it is assumed that NO BLANKS c - come through in card...only numbers (0-9) c - and maybe ONE decimal, and maybe ONE c - leading negative sign(in the leftmost c - character). ccccccccccccccccccccccccccccccccccccccccccccccccccc character*(*) card integer*4 it,iv real*8 xv,xsum il = len(card) iv = -999 xv = -999.d0 idec = -1 c - Determine if the number is negative if(card(1:1).eq.'-')then iscale = -1 ifirst = 2 else iscale = +1 ifirst = 1 endif c - Determine if there is a decimal point, and if c - so, where it falls in the character. do 1 i=ifirst,il if(card(i:i).eq.'.')then if(idec.eq.-1)then idec = i else stop 'bad value in val' endif endif 1 continue c - Interpret the integer, and make pos/neg as needed if(idec.eq.-1)then isum = 0 it = 0 do 2 i=ifirst,il iexp = (il-i) read(card(i:i),'(i1)')idum isum = isum + idum * 10**iexp 2 continue iv = isum iv = iv * iscale c - Interpret the real, and make pos/neg as needed else xsum = 0.d0 it = 1 c - Real left of decimal do 3 i=ifirst,idec-1 iexp = (idec-1)-i read(card(i:i),'(i1)')idum xsum = xsum + idum * 10**iexp 3 continue c - Real right of decimal do 4 i=idec+1,il iexp = (idec )-i read(card(i:i),'(i1)')idum iiexp = -iexp xsum = xsum + dble(idum) / 10**(iiexp) 4 continue xv = xsum xv = xv * iscale endif return end integer*4 function irev4(k) c c reverse the byte order of a four byte integer c thus changing big endian representation to little endian c (and vice-versa) c parameter (dlen=4) byte af1(dlen),af2(dlen) integer*4 k,j1,j2 equivalence (af1(1),j1),(af2(1),j2) c j1=k do 10 i=1,dlen 10 af2(i)=af1(dlen-i+1) irev4=j2 return end real*4 function frev4(f) c c reverse the byte order of a four byte float c thus changing big endian representation to little endian c (and vice-versa) c parameter (dlen=4) byte af1(dlen),af2(dlen) real*4 f,f1,f2 equivalence (af1(1),f1),(af2(1),f2) c f1=f do 10 i=1,dlen 10 af2(i)=af1(dlen-i+1) frev4=f2 return end real*8 function frev8(f) c c reverse the byte order of a eight byte float c thus changing big endian representation to little endian c (and vice-versa) c parameter (dlen=8) byte af1(dlen),af2(dlen) real*8 f,f1,f2 equivalence (af1(1),f1),(af2(1),f2) c f1=f do 10 i=1,dlen 10 af2(i)=af1(dlen-i+1) frev8=f2 return end