c - Program intd, Version 1.3, November 2, 1999 c - "Interpolate from Deflection files" c - Program will do the following: c - 1) Take two or more associated deflcetion files c - in binary (*.bin) format and compute, via c - biquadratic interpolation, the c - N/S and E/W components of the deflection of the c - vertical at any given c - latitude/longitude combination. In addition, c - the horizontal Laplace azimuth correction is computed. c - This interpolation can be done for c - one or more points, interactively c - or through input/output files. c - Note that PAIRS of files must exist together for this c - program to work (i.e. For each Xi grid, there must c - be an associated Eta grid) c - Supported Deflection models as of version 1.3: c 1) DEFLEC99 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 24, 1999 Dru A. Smith c 1.1 Oct 7, 1999 Dru A. Smith c 1.2 Oct 12, 1999 Dru A. Smith c 1.3 Nov 2, 1999 Dru A. Smith 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 86 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) Re-fixed #1 above c 2) Fixed a bug that was putting negative signs into the Blue c Book output file c c - Modifications for version 1.3: 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 c For further information, questions, or comments: c Dru A. Smith, Ph.D. c NOAA, National Geodetic Survey, N/NGS5 c U.S.A. c Phone : 301-713-3202 c Fax : 301-713-4172 c e-mail : dru@ngs.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*14 pcdirex c - Keep the Xi and Eta files separate but parallel, c - as far as the index goes integer*4 line(25),linx(25) character*256 fnamx(25),fname(25) integer*4 iosx(25),iose(25) c - No need for a duplication of these arrays...on the c - fly just make sure each header of Xi agrees with c - each header of Eta real*8 glamn(25),glamx(25),glomn(25),glomx(25) real*8 dla(25),dlo(25) integer*4 nla(25),nlo(25) integer*4 ikind(25) c - Dummy header variables for on-the-fly check real*8 glamnd,glomnd,dlad,dlod integer*4 nlad,nlod,ikindd real*4 valx,vale,valh character*80 card,card2,card85 character*40 b40 character*33 b33 character*23 sname character*17 b17 character*40 xlatc,xlonc character*1 dirxi,direta real*8 xlat,xlon real*8 pi logical ne,se,we,ee c - Variables for statistics of the run integer*4 kount,keep real*8 xn,xs,xw,xe c - Xi stats real*8 maxx,minx,maxxlat,maxxlon,minxlat,minxlon real*8 avex,stdx,rmsx c - Eta stats real*8 maxe,mine,maxelat,maxelon,minelat,minelon real*8 avee,stde,rmse c - Horizontal Laplace correction stats real*8 maxh,minh,maxhlat,maxhlon,minhlat,minhlon real*8 aveh,stdh,rmsh 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/ iosx,iose common /edges/ ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Zero out the statistics for this run ccccccccccccccccccccccccccccccccccccccccccccccccccc keep = 0 kount = 0 xn = -99999.d0 xs = +99999.d0 xw = 99999.d0 xe = -99999.d0 avex = 0.d0 stdx = 0.d0 rmsx = 0.d0 maxx = -99999999.d0 minx = +99999999.d0 avee = 0.d0 stde = 0.d0 rmse = 0.d0 maxe = -99999999.d0 mine = +99999999.d0 aveh = 0.d0 stdh = 0.d0 rmsh = 0.d0 maxh = -99999999.d0 minh = +99999999.d0 b40 = ' '// * ' ' b33 = ' ' b17 = ' ' pi = 2.d0*dasin(1.d0) 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/deflection values c - linx = set of Xi files (*.bin) c - line = set of Eta files (*.bin) ccccccccccccccccccccccccccccccccccccccccccccccccccc do 1 i=1,25 linx(i) = 10 + i line(i) = 35 + i 1 continue linn = 1 lout = 10 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Which deflection model? ccccccccccccccccccccccccccccccccccccccccccccccccccc 108 write(6,101) 101 format(/,1x,70('-'),/, * ' Which deflection model do you wish to use?',//, * ' 1 = DEFLEC99 ',/, * ' ',/, * ' 99 = END PROGRAM',//, * ' -> ',$) read(5,*)imodel if(imodel.eq.99)stop if(imodel.lt.1 .or. imodel.gt.1)goto 108 ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Which directory are the deflection files in? ccccccccccccccccccccccccccccccccccccccccccccccccccc pcdirex = char(39)//'C:'//char(92)//'DEFLEC99'//char(92)//char(39) write(6,102)pcdirex 102 format(/,1x,70('-'),/, * ' What is the **FULL** directory name (including', * ' trailing slashes) ',/, * ' where the deflection (*.bin) files may be', * ' found?',/, * 5x,' (Unix Example: ''/export/home/deflec99/'') ',/, * 5x,' (PC Example: ',a14,') ',/, * 5x,' Hit to default to this directory',//, * ' -> ',$) 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 (iosx/iose) c - of which files are open and which are not. c - Abort inside this subroutine if companion c - files are not found for any file. ccccccccccccccccccccccccccccccccccccccccccccccccccc call files(imodel,dirnam,nfiles,fnamx,fname,linx,line,nff) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Read the headers of all deflection files which c - were opened, and store that information. c - Compute the max lat/lon from the header information. c - Do on-the-fly checks of Xi headers against c - Eta headers. ccccccccccccccccccccccccccccccccccccccccccccccccccc do 115 i=1,nfiles / 2 if(iosx(i).eq.0)then c - Xi read(linx(i),rec= 1)glamnx(1) read(linx(i),rec= 2)glamnx(2) read(linx(i),rec= 3)glomnx(1) read(linx(i),rec= 4)glomnx(2) read(linx(i),rec= 5)dlax(1) read(linx(i),rec= 6)dlax(2) read(linx(i),rec= 7)dlox(1) read(linx(i),rec= 8)dlox(2) read(linx(i),rec= 9)nla(i) read(linx(i),rec=10)nlo(i) read(linx(i),rec=11)ikind(i) glamn(i) = glamny glomn(i) = glomny dla(i) = dlay dlo(i) = dloy c - Eta read(line(i),rec= 1)glamnx(1) read(line(i),rec= 2)glamnx(2) read(line(i),rec= 3)glomnx(1) read(line(i),rec= 4)glomnx(2) read(line(i),rec= 5)dlax(1) read(line(i),rec= 6)dlax(2) read(line(i),rec= 7)dlox(1) read(line(i),rec= 8)dlox(2) read(line(i),rec= 9)nlad read(line(i),rec=10)nlod read(line(i),rec=11)ikindd glamnd = glamny glomnd = glomny dlad = dlay dlod = dloy c - On the fly check: if(glamnd .ne. glamn(i))stop 70001 if(glomnd .ne. glomn(i))stop 70002 if(dlad .ne. dla(i) )stop 70003 if(dlod .ne. dlo(i) )stop 70004 if(nlad .ne. nla(i) )stop 70005 if(nlod .ne. nlo(i) )stop 70006 if(ikindd .ne. ikind(i))stop 70007 c - Compute lat max and lon max 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 Deflections) Type 1',/, * ' 2 = Free Format (For Deflections) Type 2',/, * ' 3 = NGS Blue Book Format(*80* and *85* 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 DEFLECTIONS) 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 - deflection values to -999 and skip the interpolation if(xlat.eq.-999.d0 .or. xlon.eq.-999.d0)then valx = -999.d0 vale = -999.d0 valh = -999.d0 goto 130 endif c - Force longitude to be positive East if(.not.poseast)xlon = 360.d0 - xlon c - Find which deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = 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, * valx,vale,valh) c - Go back and get another value goto 116 c - When finished, give a little report to the c - screen and end program. 114 avex = avex/keep rmsx = sqrt(rmsx/keep) avee = avee/keep rmse = sqrt(rmse/keep) aveh = aveh/keep rmsh = sqrt(rmsh/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) stdx = sqrt(fact*abs(rmsx**2 - avex**2)) stde = sqrt(fact*abs(rmse**2 - avee**2)) stdh = sqrt(fact*abs(rmsh**2 - aveh**2)) else stdx = 0 stde = 0 stdh = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * avex,stdx,minx,minxlat,minxlon,maxx,maxxlat,maxxlon, * avee,stde,mine,minelat,minelon,maxe,maxelat,maxelon, * aveh,stdh,minh,minhlat,minhlon,maxh,maxhlat,maxhlon else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * avex,stdx,minx,minxlat,360.d0-minxlon,maxx, * maxxlat,360.d0-maxxlon, * avee,stde,mine,minelat,360.d0-minelon,maxe, * maxelat,360.d0-maxelon, * aveh,stdh,minh,minhlat,360.d0-minhlon,maxh, * maxhlat,360.d0-maxhlon endif stop 201 format(/,1x,70('-'),/, * ' FINAL REPORT: ',/, * ' Points Input / Points Kept : ',2x,i8,2x,i8,/, * ' North/South/West/East Bounds: ',2x,4(f10.6,1x),//, * ' Basic Statistics: ',/, * ' Val AVE Sigma Minim. Lat of Min Lon ', * 'of Min Maxim. Lat of Max Lon of Max',/, * ' (sec) (sec) (sec) (dec deg) ', * '(dec deg) (sec) (dec deg) (dec deg)',/, * ' Xi ',f7.2,1x,f6.2,1x,f7.2,1x,f10.6,1x, * f10.6,1x,f7.2,1x,f10.6,1x,f10.6,/, * ' Eta ',f7.2,1x,f6.2,1x,f7.2,1x,f10.6,1x, * f10.6,1x,f7.2,1x,f10.6,1x,f10.6,/, * ' Lap ',f7.2,1x,f6.2,1x,f7.2,1x,f10.6,1x, * f10.6,1x,f7.2,1x,f10.6,1x,f10.6,//) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - INPUT FILE, FREE FORMAT (FOR DEFLECTIONS) 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 - deflection values to -999 and skip the interpolation if(xlat.eq.-999.d0 .or. xlon.eq.-999.d0)then valx = -999.d0 vale = -999.d0 valh = -999.d0 goto 330 endif c - Force longitude to be positive East if(.not.poseast)xlon = 360.d0 - xlon c - Find which deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = 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, * valx,vale,valh) c - Go back and get another value goto 119 c - When finished, give a little report to the c - screen and end program. 118 avex = avex/keep rmsx = sqrt(rmsx/keep) avee = avee/keep rmse = sqrt(rmse/keep) aveh = aveh/keep rmsh = sqrt(rmsh/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) stdx = sqrt(fact*abs(rmsx**2 - avex**2)) stde = sqrt(fact*abs(rmse**2 - avee**2)) stdh = sqrt(fact*abs(rmsh**2 - aveh**2)) else stdx = 0 stde = 0 stdh = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * avex,stdx,minx,minxlat,minxlon,maxx,maxxlat,maxxlon, * avee,stde,mine,minelat,minelon,maxe,maxelat,maxelon, * aveh,stdh,minh,minhlat,minhlon,maxh,maxhlat,maxhlon else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * avex,stdx,minx,minxlat,360.d0-minxlon,maxx, * maxxlat,360.d0-maxxlon, * avee,stde,mine,minelat,360.d0-minelon,maxe, * maxelat,360.d0-maxelon, * aveh,stdh,minh,minhlat,360.d0-minhlon,maxh, * maxhlat,360.d0-maxhlon 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 *85* c - record. If so, overwrite the deflection fields in c - it with interpolated values. If not, c - create an *85* record. c - Read one line and 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 *85* 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 *85* 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 85 record are 81, and 82 records (83 and 84 c - records will be deleted). Also, 86 records must c - be allowed to follow the 85 record 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*, or *82* 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')then if(outfil)then write(lout,136)card2 endif goto 140 endif c - POSSIBILITY 3: We find an *85* record c - ACTION: a) If it is the wrong *85 record, delete c it and make a NEW *85* record and then c go back and get another "card2" value c - b) Modify it if it's the right *85* record if(card2(8:9).eq.'85')then if(card2(11:14).eq.card(11:14))then c - Arriving here, we've found the *85* record associated with c - the *80* record c - Find which deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = xlon endif endif c - Modify the *85* record and write it out c - Find the right deflection code c - DEFLEC99 if(imodel.eq.1)then card2(21:29) = ' DEFLEC99' card2(30:62) = b33 endif ixi = nint(valx*100) ieta = nint(vale*100) if(outfil)then if(ixi.ge.0)dirxi='N' if(ixi.lt.0)dirxi='S' if(ieta.ge.0)direta='E' if(ieta.lt.0)direta='W' write(lout,936)card2(1:62),abs(ixi),dirxi,'000', * abs(ieta),direta,'000' endif c - And now go back and look for the next "card" value c - (Any *86* records which follow this *85* record c - will be passed into the output file) goto 137 c - If this is NOT the associated *85* record, delete c - it and make a new *85* record else c - Find which deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = xlon endif endif c - Make a NEW *85* record and write it out c - Find the right deflection code c - DEFLEC99 card85 = b40//b40 if(imodel.eq.1)then card85(21:29) = ' DEFLEC99' card85(30:62) = b33 endif ixi = nint(valx*100) ieta = nint(vale*100) if(outfil)then if(ixi.ge.0)dirxi='N' if(ixi.lt.0)dirxi='S' if(ieta.ge.0)direta='E' if(ieta.lt.0)direta='W' write(lout,936)card85(1:62),abs(ixi),dirxi,'000', * abs(ieta),direta,'000' endif c - And now go back and look for the next "card" value c - (Any *86* records which follow this *85* record c - will be passed into the output file) goto 137 endif endif c - POSSIBILITY 4: We find something besides an 81,82, 83,84 c or 85 record c - ACTION: Create a new *85* 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' )then c - Find which deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = xlon endif endif c - Make a NEW *85* record ixi = nint(valx*100) ieta = nint(vale*100) card85 = b40//b40 card85(1:6) = card(1:6) card85(7:10) = '*85*' card85(11:14) = card(11:14) c - Find the right deflection comments c - DEFLEC99 if(imodel.eq.1)then card85(21:29) = ' DEFLEC99' card85(30:62) = b33 endif c - ...and write out the *85* record if(outfil)then if(ixi.ge.0)dirxi='N' if(ixi.lt.0)dirxi='S' if(ieta.ge.0)direta='E' if(ieta.lt.0)direta='W' write(lout,936)card85(1:62),abs(ixi),dirxi,'000', * abs(ieta),direta,'000' endif 936 format(a62,i5,a1,a3,i5,a1,a3) c - Now put the latest 'card2' record into 'card' and go back to c - the top to search...if 'card2' is an *86* record, then c - it will be passed unchanged into the output file. card = card2 goto 153 endif c - POSSIBILITY 5: We find the EOF c - ACTION: Create a new *85* based on the previous *80* record c - and then go to the 'end report' phase c - Find which deflection 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 valx = -999 vale = -999 valh = -999 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = xlon endif endif c - Make a NEW *85* record ixi = nint(valx*100) ieta = nint(vale*100) card85 = b40//b40 card85(1:6) = card(1:6) card85(7:10) = '*85*' card85(11:14) = card(11:14) c - Find the right deflection comments c - DEFLEC99 if(imodel.eq.1)then card85(21:29) = ' DEFLEC99' card85(30:62) = b33 endif c - ...and write out the *85* record if(outfil)then if(ixi.ge.0)dirxi='N' if(ixi.lt.0)dirxi='S' if(ieta.ge.0)direta='E' if(ieta.lt.0)direta='W' write(lout,936)card85(1:62),abs(ixi),dirxi,'000', * abs(ieta),direta,'000' 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 avex = avex/keep rmsx = sqrt(rmsx/keep) avee = avee/keep rmse = sqrt(rmse/keep) aveh = aveh/keep rmsh = sqrt(rmsh/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) stdx = sqrt(fact*abs(rmsx**2 - avex**2)) stde = sqrt(fact*abs(rmse**2 - avee**2)) stdh = sqrt(fact*abs(rmsh**2 - aveh**2)) else stdx = 0 stde = 0 stdh = 0 endif write(6,201)kount,keep,xn,xs,xw,xe, * avex,stdx,minx,minxlat,minxlon,maxx,maxxlat,maxxlon, * avee,stde,mine,minelat,minelon,maxe,maxelat,maxelon, * aveh,stdh,minh,minhlat,minhlon,maxh,maxhlat,maxhlon 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 deflection files 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 valx = -999.d0 vale = -999.d0 valh = -999.d0 else c - Otherwise, do the interpolation c - Xi first call interp(xlat,xlon,linx,k,valx) c - Eta next call interp(xlat,xlon,line,k,vale) c - Finally the Horizontal Laplace correction valh = -tan(xlat*pi/180.d0) * vale keep = keep + 1 if(xlat.gt.xn)xn = xlat if(xlat.lt.xs)xs = xlat if(xlon.gt.xe)xe = xlon if(xlon.lt.xw)xw = xlon avex = avex + valx rmsx = rmsx + valx*valx if(valx.lt.minx)then minx = valx minxlat = xlat minxlon = xlon endif if(valx.gt.maxx)then maxx = valx maxxlat = xlat maxxlon = xlon endif avee = avee + vale rmse = rmse + vale*vale if(vale.lt.mine)then mine = vale minelat = xlat minelon = xlon endif if(vale.gt.maxe)then maxe = vale maxelat = xlat maxelon = xlon endif aveh = aveh + valh rmsh = rmsh + valh*valh if(valh.lt.minh)then minh = valh minhlat = xlat minhlon = xlon endif if(valh.gt.maxh)then maxh = valh maxhlat = xlat maxhlon = xlon endif endif c - Finally, write out the record to screen and possibly c - to an output file. Use Free Format (For Deflections) Type 1 c - output card = sname//b17//b40 call ff1out(iinput,outfil,lout,card,poseast,xlat,xlon, * valx,vale,valh) 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 avex = avex/keep rmsx = sqrt(rmsx/keep) avee = avee/keep rmse = sqrt(rmse/keep) aveh = aveh/keep rmsh = sqrt(rmsh/keep) if(keep.gt.1)then fact = dble(keep)/dble(keep-1) stdx = sqrt(fact*abs(rmsx**2 - avex**2)) stde = sqrt(fact*abs(rmse**2 - avee**2)) stdh = sqrt(fact*abs(rmsh**2 - aveh**2)) else stdx = 0 stde = 0 stdh = 0 endif if(poseast)then write(6,201)kount,keep,xn,xs,xw,xe, * avex,stdx,minx,minxlat,minxlon,maxx,maxxlat,maxxlon, * avee,stde,mine,minelat,minelon,maxe,maxelat,maxelon, * aveh,stdh,minh,minhlat,minhlon,maxh,maxhlat,maxhlon else write(6,201)kount,keep,xn,xs,360.d0-xw,360.d0-xe, * avex,stdx,minx,minxlat,360.d0-minxlon,maxx, * maxxlat,360.d0-maxxlon, * avee,stde,mine,minelat,360.d0-minelon,maxe, * maxelat,360.d0-maxelon, * aveh,stdh,minh,minhlat,360.d0-minhlon,maxh, * maxhlat,360.d0-maxhlon endif stop endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - End Of Main Program ccccccccccccccccccccccccccccccccccccccccccccccccccc end c c --------------------------------------------------------------- c subroutine files(imodel,dirnam,nfiles,fnamx,fname,linx,line,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 linx(25),line(25) character*256 fnamx(25),fname(25) integer*4 iosx(25),iose(25) 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/ iosx,iose 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 DEFLEC99 file names cccccccccccccccccccccccccccccccccccccc 2 if(imodel.eq.1)then nfiles = 28 c - First 8 areas are CONUS do 101 i=1,8 write(cval,109)i fnamx(i) = dirnam(1:ilen)//'x1999u'//cval//suf fname(i) = dirnam(1:ilen)//'e1999u'//cval//suf 101 continue c - Next 1 file is HAWAII do 102 i=9,9 write(cval,109)i-8 fnamx(i) = dirnam(1:ilen)//'x1999h'//cval//suf fname(i) = dirnam(1:ilen)//'e1999h'//cval//suf 102 continue c - Next 1 file is PR/VI do 103 i=10,10 write(cval,109)i-9 fnamx(i) = dirnam(1:ilen)//'x1999p'//cval//suf fname(i) = dirnam(1:ilen)//'e1999p'//cval//suf 103 continue c - Next 4 files are Alaska do 104 i=11,14 write(cval,109)i-10 fnamx(i) = dirnam(1:ilen)//'x1999a'//cval//suf fname(i) = dirnam(1:ilen)//'e1999a'//cval//suf 104 continue endif ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Now open all the files ccccccccccccccccccccccccccccccccccccccccccccccccccc do 201 i=1,nfiles / 2 open(linx(i),file=fnamx(i),status='old',form='unformatted', * access='direct',recl=4,iostat=iosx(i)) open(line(i),file=fname(i),status='old',form='unformatted', * access='direct',recl=4,iostat=iose(i)) if(iosx(i).eq.0)then write(6,203) write(6,204) fnamx(i)(1:ilen+12) endif if(iose(i).eq.0)then write(6,203) write(6,204) fname(i)(1:ilen+12) endif c - Check for companion files if(iosx(i).eq.0 .and. iose(i).ne.0)then write(6,207) write(6,204)fnamx(i)(1:ilen+12) stop endif if(iose(i).eq.0 .and. iosx(i).ne.0)then write(6,207) write(6,204)fname(i)(1:ilen+12) stop endif 201 continue 203 format(' *** Opening File: ',$) 204 format(a) 207 format(' *** ERROR(207): Can not find companion to file: ',$) 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 deflection files. c - No need to check both, as the check for c - companion files was done already ccccccccccccccccccccccccccccccccccccccccccccccccccc nff = 0 do 202 i=1,nfiles / 2 if(iosx(i).eq.0)nff = nff + 2 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 Deflections) Type 1',/, * ' 2 = Free Format (For Deflections) 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 Deflections) Type 1 ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iformx.eq.1)then write(6,3) 3 format(/,1x,70('-'),/, * ' 1 = Free Format (For Deflections) Type 1 : ',/, * ' This is an ASCII format, where the first 23 characters',/, * ' of each record may contain the station name or be blank.',/, * ' The last 40 characters (41-80) of the record 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,12x,'AIRPORT SC1',17x,' 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 Deflections) Type 2 ccccccccccccccccccccccccccccccccccccccccccccccccccc if(iformx.eq.2)then write(6,5) 5 format(/,1x,70('-'),/, * ' 2 = Free Format (For Deflections) Type 2 : ',/, * ' This is an ASCII format, where the last 23 characters', * ' (58-80)',/, * ' 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,25x,'BM1027',/, * ' 45.3 282',15x,25x,'XXX8063',/, * ' 47 18.55 252 14',15x,25x,8x,'AIRPORT SC1',/, * ' 30 17 270 59',/, * ' 48 14 17 239 52 18.753',10x,25x,7x,'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 *85* (Deflection)', /, + ' Records. INTD uses information from the', + ' *80* and *85* 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', /, + ' *85* record is modified by INTD - the *80*', + ' records and all other records are') WRITE (6, 511) 511 FORMAT (' passed through without change to the output file.', + ' Information in columns', + ' 21-80 in the *85* record is', /, + ' created by INTD.',/) 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 *85* records', + ' from a Blue Book', /, + ' file:', //, + '004560*80*0096KNOXVILLE CT HSE ', + '411906578 N0930548534 W 277 MIA33', /, + '004565*85*0096 any deflections 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,' INTD PROGRAM',/, * 10x,' (Interpolate from Deflection files). ',/) write(6,2) 2 format( * 10x,' VERSION 1.3',/, * 10x,' November 2, 1999',/, * 10x,' Dru A. Smith, Ph.D.',//, * 10x,' (Hit RETURN to continue)') read(5,'(a)')keyb ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Disclaimer ccccccccccccccccccccccccccccccccccccccccccccccccccc WRITE (6,932) 932 format(/,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, (For c - Deflections) Type 1" record (card) and extracts c - the latitude and longitude (xlat,xlon) in c - 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 decimals '.', 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 - deflection files will be used to interpolate c - to a point at xlat/xlon. The returned c - value of "k" means the "kth" pair of files will c - be used. c - c - For the DEFLEC99 models: c Alaska and CONUS overlap, and this code c - forces a "CONUS wins" scenario when interpolating c - the deflections at points in the overlap region. ccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 xlat,xlon integer*4 nfiles,k,imodel real*8 glamn(25),glamx(25),glomn(25),glomx(25) real*8 dla(25),dlo(25) integer*4 nla(25),nlo(25),ikind(25) integer*4 rank(25) integer*4 iosx(25),iose(25) 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/ iosx,iose common /edges/ ne,se,we,ee ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Spin through all *open* pairs of files, and c - *RANK* them...the pair of files 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 pair's area c - 1 = Point lies in this pair, but at a corner c - 2 = Point lies in this pair, but at an edge c - 3 = Point lies in this pair, 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 / 2 if(iosx(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 pair of grids 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 DEFLEC99 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. iosx(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 pair, 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 pair of files which c - had a rank of 3. So now, see if we have any rank 2 c - or rank 1 pairs to use. do 101 i=1,nfiles / 2 if(iosx(i).ne.0)goto 101 if(rank(i).eq.2)then k = i return endif 101 continue do 102 i=1,nfiles / 2 if(iosx(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, * valx,vale,valh) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to write out a "Free Format, (For c - Deflections) Type 1" record ccccccccccccccccccccccccccccccccccccccccccccccccccc logical outfil logical poseast integer*4 lout,iinput character*80 card character*40 fcard real*8 xlat,xlon real*4 valx,vale,valh integer*4 deglat,minlat,deglon,minlon real*4 seclat,seclon fcard = card(1:23) 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,valx,vale,valh endif if(iinput.eq.1)write(6,98)fcard,deglat,minlat,seclat, * deglon,minlon,seclon,valx,vale,valh 98 format(/,' Your Result:',/, * 1x,23x,1x,' latitude ',1x, * ' longitude ',1x,' Xi ',1x,' Eta ',1x, * 'Hor Lap',/, * ' Station Name',9x,1x,'ddd mm ss.sssss',1x, * 'ddd mm ss.sssss'1x,'arc-sec',1x,'arc-sec', * 1x,'arc-sec',/, * 1x,a23,1x,i3,1x,i2,1x,f8.5, * 1x,i3,1x,i2,1x,f8.5, * 1x,f7.2,1x,f7.2,1x,f7.2,/) 99 format(a23,1x,i3,1x,i2,1x,f8.5, * 1x,i3,1x,i2,1x,f8.5, * 1x,f7.2,1x,f7.2,1x,f7.2) return end c c --------------------------------------------------------------- c subroutine ff2out(outfil,lout,card,poseast,xlat,xlon, * valx,vale,valh) ccccccccccccccccccccccccccccccccccccccccccccccccccc c - Subroutine to write out a "Free Format, (For c - Deflections) Type 2" record ccccccccccccccccccccccccccccccccccccccccccccccccccc logical outfil logical poseast integer*4 lout character*80 card character*23 ecard real*8 xlat,xlon real*4 valx,vale,valh integer*4 deglat,minlat,deglon,minlon real*4 seclat,seclon ecard = card(58: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, * valx,vale,valh,ecard endif 99 format(1x,i3,1x,i2,1x,f8.5, * 1x,i3,1x,i2,1x,f8.5, * 1x,f7.2,1x,f7.2,1x,f7.2,a23) 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(25) integer*4 k real*4 val real*8 glamn(25),glamx(25),glomn(25),glomx(25) real*8 dla(25),dlo(25) integer*4 nla(25),nlo(25) integer*4 ikind(25) real*4 f1,f2,f3,f4,f5,f6,f7,f8,f9 logical ne,se,we,ee 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 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 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