program deflec *********************************************************************** * * * program : deflec * * * * purpose: computation program to interpolate deflection * * values from a grid of regularly-spaced estimates. * * the input requires latitude and longitude values * * referrenced to the north american datum of 1983 (nad83),* * for the conterminous U.S. * * input data referrenced to the north american datum of * * 1927 (nad27) should be converted to nad 83 via nadcon * * prior to using this program, for the USA. NAD 27 is * * acceptable for MEXICO and the CARIBBEAN. * * * * the actual computation is performed as an interpolation * * from a regularly-spaced grid of points obtained from * * the deflec96 model computed by dr. dennis g. milbert * * and dr. dru a. smith. * * * * VERSION 3.10: * * MODELS DMEX97 (computed by Dr. Dru A. Smith) and * * DCAR97 (computed by Dru A. Smith and * * Howard J. Small (NIMA)) have been added * * NOTE THAT THE MODELS DMEX97 AND DCAR97 ARE * * ALL RELATIVE TO THE GRS-80 ELLIPSOID, CENTERED * * AT THE ORIGIN OF THE ITRF94(1996.0) REFERENCE FRAME. * * YOU CAN *STILL* INPUT LAT/LON IN NAD83 AND EXPECT TO * * GET THE SAME VALUE AS IF YOU HAD USED GRS-80/ITRF94 * * LAT/LON VALUES, HOWEVER IT *MUST* BE POINTED OUT THAT * * THE GEOID HEIGHTS OF THOSE MODELS ARE *NOT* GEOID * * HEIGHTS ABOVE NAD83! (THE REASON THIS WORKS IS BECAUSE * * THE HORIZONTAL SHIFT FROM NAD83 TO ITRF94 IS ABOUT 1 OR * * 2 METERS, WHICH WILL HAVE *NO* IMPACT ON THE GEOID * * HEIGHT VALUE; HOWEVER THE VERTICAL SHIFT FROM ITRF94 * * TO NAD83 IS ALSO 1-2 METERS, SO IT IS IMPORTANT TO * * REMEMBER WHICH ELLIPSOID THE GEOIDS REFER TO. * * * * * * the interpolation is accomplished by a locally fit * * biquadratic function. * * * * the polynomial surface is fit to nine data points * * defining the area surrounding the (x,y) pair * * where the interpolation is to take place. * * * * the program requires that the user specify: * * * * 1) the name of an output file * * * * 2) the name of an input file (if available). * * * * * * this program allows for either ngs standard horizontal * * data formats as specified in the fgcc publication, * * commonly known as the 'horizontal blue book' (see * * subroutines type3, type4, type5), or in a generic file * * format (see subroutine type1 or subroutine type2). * * * * the code can be easily modified to accommodate custom * * file specifications by modifying subroutines: endrep, * * getpt, iparms, wrtpt, and (optionally) fhelp. * * * * * * version code: 3.00 (20-oct-96) * * version code: 3.10 ( 8-apr-97) * * * * version date: may 7, 1997 * * july 16, 1997 -- retained version # as 3.10 * * * * added support for inexact header in GRIDS() * * increaded MXAREA to 18 * * * * authors: dru a. smith, ph.d. * * dennis g. milbert, ph.d. * * systems development division * * national geodetic survey, nos, noaa * * silver spring, md 20910 * * * * original authors: (software) * * author: warren t. dewhurst, ph.d. (software) * * lieutenant commander, noaa * * alice r. drew (software) * * senior geodesist, horizontal network branch * *********************************************************************** *********************************************************************** * * * 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. * * * *********************************************************************** * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) real vrsion integer mxarea parameter (vrsion = 3.10e0, mxarea =18) double precision xsmall, xbig, ysmall, ybig double precision axii, vxii double precision sdxii, sdxii2 real xii, smxii, bgxii double precision aeta, veta double precision sdeta, sdeta2 real eta, smeta, bgeta integer ipage, itype, nconv, ifile character*80 name logical page, nodata, screen real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) * v 3.10: differentiate between deflec96, dmex97 and dcar97, das logical ldmex,ldcar common/scienc/ldmex,ldcar ********************** * initialize variables ********************** call initl (screen, page, name, ipage, itype, + smxii, bgxii, axii, vxii, sdxii, sdxii2, + smeta, bgeta, aeta, veta, sdeta, sdeta2, + xsmall, xbig, ysmall, ybig) ************************** * print header information ************************** call headr (vrsion) *********************** * open deflec data files *********************** ldmex=.false. ldcar=.false. call ngrids (nodata) if (nodata) goto 9999 ****************************************************** * request for the needed variable values from the user ****************************************************** call iparms (itype, screen) ******************************** * loop once for each computation ******************************** call mloop (nconv, ipage, itype, vrsion, + xii, sdxii, sdxii2, smxii, bgxii, + eta, sdeta, sdeta2, smeta, bgeta, + xsmall, xbig, ysmall, ybig, + page, screen, name) *********************************************** * finished with all computations - write report *********************************************** call endrep (ipage, nconv, itype, + xii, axii, vxii, sdxii, sdxii2, smxii, bgxii, + eta, aeta, veta, sdeta, sdeta2, smeta, bgeta, + xsmall, xbig, ysmall, ybig) ***************** * close all files ***************** do 1010 ifile = 1, 2*narea close ( luarea(ifile) ) 1010 continue close (nin, status='keep') close (nout, status='keep') close (napar, status='keep') 9999 stop end subroutine askpt (nconv, name, idla, imla, sla, + idlo, imlo, slo, xpt, ypt, eof, nopt) * interactively ask for the name and location of a point * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) character*40 b40 parameter (b40 = ' ') double precision xpt, ypt, rdla, rdlo, dcard real sla, slo, rmla, rmlo, rcard integer nconv integer idla, imla, idlo, imlo integer iflag1, iflag2, n1, n2 integer leng, ierr, ios character*80 name character*40 ans, dum logical eof, nopt integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) data iflag1 /1/, iflag2 /2/ name = ' ' write (luout,*) ' What is the NAME for this station or', + ' point?' read (luin,'(a80)') name ********** * latitude ********** if (nconv .eq. 0) then write (luout,110) 110 format (/, ' Latitudes and longitudes may be entered', + ' in three formats:', /, + ' (1) degrees, minutes, and decimal seconds, OR', /, + ' (2) degrees, and decimal minutes OR', /, + ' (3) decimal degrees.', /, + ' Degrees, minutes and seconds may be separated', + ' by blanks or commas.', /, + ' a latitude or longitude of 0 will end data', + ' entry.', /) WRITE (LUOUT,1111) 1111 FORMAT (' *************** Which latitude/longitude to use? ', + ' ***************',/, + ' You may enter latitudes and longitudes in either the', + ' NAD 83 or the',/, + ' ITRF94/GRS-80 system, with no change in the interpolated', + ' deflection value.',/, + ' Do not use NAD 27 in the United States.',/, + ' Transform to NAD 83 first,', + ' using NADCON software.',/) endif write (luout,*) ' What is its geodetic latitude?' write (luout,*) ' ' read (luin,170,err=9930,iostat=ios) ans 170 format (a40) if (ans .eq. b40) goto 9999 dum = ans call nblank (dum, iflag2, n2) leng = n2 rdla = dcard( dum(1:n2), leng, ierr ) if (ierr .ne. 0) goto 9950 if (leng .gt. 0) then rmla = rcard( dum, leng, ierr ) if (ierr .ne. 0) goto 9950 if (leng .gt. 0) then sla = rcard( dum, leng, ierr ) if (ierr .ne. 0) goto 9950 else sla = 0.e0 endif else rmla = 0.e0 sla = 0.e0 endif if ( (rdla .eq. 0.d0) .and. (rmla .eq. 0.e0) .and. + (sla .eq. 0.e0) ) goto 9999 * check for illogical values if (rdla .lt. 0.d0) goto 9940 if (rdla .gt. 90.d0) goto 9950 if (rmla .lt. 0.e0 .or. rmla .gt. 60.e0) goto 9950 if ( sla .lt. 0.e0 .or. sla .gt. 60.e0) goto 9950 *********** * longitude *********** write (luout,*) ' What is its geodetic longitude?', + ' (Longitude is positive WEST.)' write (luout,*) ' ' read (luin,170,err=9930,iostat=ios) ans if (ans .eq. b40) goto 9999 dum = ans call nblank (dum, iflag2, n2) leng = n2 rdlo = dcard( dum(1:n2), leng, ierr ) if (ierr .ne. 0) goto 9960 if (leng .gt. 0) then rmlo = rcard( dum, leng, ierr ) if (ierr .ne. 0) goto 9960 if (leng .gt. 0) then slo = rcard( dum, leng, ierr ) if (ierr .ne. 0) goto 9960 else slo = 0.e0 endif else rmlo = 0.e0 slo = 0.e0 endif if ( (rdlo .eq. 0.d0) .and. (rmlo .eq. 0.e0) .and. + (slo .eq. 0.e0) ) goto 9999 * check for illogical values if (rdlo .lt. 0.d0) goto 9940 if (rdlo .gt. 360.d0) goto 9960 if (rmlo .lt. 0.e0 .or. rmlo .gt. 60.e0) goto 9960 if ( slo .lt. 0.e0 .or. slo .gt. 60.e0) goto 9960 * calculate decimal degrees ypt = rdla + dble(rmla)/60.d0 + dble(sla)/3600.d0 xpt = rdlo + dble(rmlo)/60.d0 + dble(slo)/3600.d0 * get degrees, minutes, seconds call hms (ypt, idla, imla, sla) call hms (xpt, idlo, imlo, slo) 9000 return * error messages 9930 continue call nblank (ans, iflag1, n1) call nblank (ans, iflag2, n2) write (luout,9935) ans(n1:n2) 9935 format (' ERROR - in the answer:', /, + 9x, '''', a, '''', /, + ' Must enter number in prescribed format!', /) nopt = .true. goto 9000 9940 write (luout,9945) call nblank (ans, iflag1, n1) call nblank (ans, iflag2, n2) write (luout,9945) ans(n1:n2) 9945 format (' ERROR - in the answer:', /, + 9x, '''', a, '''', /, + ' Latitude and longitudes must be positive!', /, + ' Longitude is positive west.', /) nopt = .true. goto 9000 9950 continue call nblank (ans, iflag1, n1) call nblank (ans, iflag2, n2) write (luout,9955) ans(n1:n2) 9955 format (' ERROR - Illogical value for latitude in the answer:', /, + ' ''', a, '''', /, + ' Latitude must be between 0 and 90 degrees.', /, + ' Minutes and seconds must be between 0', + ' and 60.', /) nopt = .true. goto 9000 9960 continue call nblank (ans, iflag1, n1) call nblank (ans, iflag2, n2) write (luout,9965) ans(n1:n2) 9965 format (' ERROR - Illogical value for longitude in the answer:',/, + ' ''', a, '''', /, + ' Longitude must be between 0 and 360 degrees.',/, + ' Minutes and seconds must be between 0', + ' and 60.', /) nopt = .true. goto 9000 9999 continue eof = .true. goto 9000 end subroutine diagrm (lu, nconv, xsmall, xbig, ysmall, ybig) * this subroutine prints out a small diagram showing the area * that was computed. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) double precision xtemp, xsmall, xbig, ysmall, ybig real slomin, slomax, slamin, slamax integer lodmin, lommin, lodmax, lommax integer ladmin, lammin, ladmax, lammax integer lu, nconv xtemp = xsmall xsmall = xbig xbig = xtemp call hms (xbig, lodmin, lommin, slomin) call hms (xsmall, lodmax, lommax, slomax) call hms (ysmall, ladmin, lammin, slamin) call hms (ybig, ladmax, lammax, slamax) write (lu,90) nconv 90 format (//, ' The total number of deflection estimates: ', i8) write (lu,100) 100 format (//, 30x, 'Region of estimation') write (lu,1000) lodmax, lommax, slomax, lodmin, lommin, slomin, + ladmax, lammax, slamax, ladmax, lammax, slamax, + ladmin, lammin, slamin, ladmin, lammin, slamin, + lodmax, lommax, slomax, lodmin, lommin, slomin 1000 format (3(/), t4, ' ', /, + t4, 'Longitude:', 8x, i4, 1x, i2.2, 1x, f6.3, 13x, i4, 1x, + i2.2, 1x, f6.3, /, + t4, 'Latitude:', 9x, i4, 1x, i2.2, 1x, f6.3, + ' ************', i4, 1x, i2.2, 1x, f6.3, /, + 5(t27, 3x, '*', t57, '*', /), + t27, 3x, '*', 9x, ' ', t57, '*', /, + t27, 3x, '*', 9x, ' data ', t57, '*', /, + t27, 3x, '*', 9x, ' points', t57, '*', /, + 5(t27, 3x, '*', t57, '*', /), + t4, 'Latitude:', 9x, i4, 1x, i2.2, 1x, f6.3, + ' ************', i4, 1x, i2.2, 1x, f6.3, /, + t4, 'Longitude:', 8x, i4, 1x, i2.2, 1x, f6.3, 13x, i4, 1x, + i2.2, 1x, f6.3, //) return end subroutine endrep (ipage, nconv, itype, + xii, axii, vxii, sdxii, sdxii2, smxii, bgxii, + eta, aeta, veta, sdeta, sdeta2, smeta, bgeta, + xsmall, xbig, ysmall, ybig) *** gather statistics and write the end-of-program report * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) double precision xsmall, xbig, ysmall, ybig double precision axii, vxii double precision sdxii, sdxii2 real xii, smxii, bgxii double precision aeta, veta double precision sdeta, sdeta2 real eta, smeta, bgeta integer ipage, nconv, itype integer lu, i logical page real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) page = .true. ipage = ipage + 1 ******************* * do the statistics ******************* if (nconv .ge. 2) then * calculate mean, variance, standard deviation for deflection * estimates (arc sec). ******************* * xi deflection (arc sec) ******************* axii = sdxii/dble(nconv) vxii = sdxii2/dble(nconv-1) - sdxii**2/dble( nconv*(nconv-1 )) if (vxii .gt. 1.0d-6) then sdxii = dsqrt(vxii) else vxii = 0.0d0 sdxii = 0.0d0 endif ******************* * eta deflection (arc sec) ******************* aeta = sdeta/dble(nconv) veta = sdeta2/dble(nconv-1) - sdeta**2/dble( nconv*(nconv-1 )) if (veta .gt. 1.0d-6) then sdeta = dsqrt(veta) else veta = 0.0d0 sdeta = 0.0d0 endif endif if (nconv .lt. 2) then axii = xii vxii = 0.0d0 sdxii = 0.0d0 aeta = eta veta = 0.0d0 sdeta = 0.0d0 endif *************************************** * print out the statistics for this job *************************************** lu = luout if (nconv .gt. 0) then ************************************************* * first report the final statistics to the screen ************************************************* call report (lu, smxii, bgxii, axii, vxii, sdxii, + smeta, bgeta, aeta, veta, sdeta, ipage, page) call diagrm (lu, nconv, xsmall, xbig, ysmall, ybig) **************************************************** * now report the final statistics to the output file **************************************************** if (itype .eq. 0) then ************************************** * interactive use only - no input file ************************************** lu = nout call report (lu, smxii, bgxii, axii, vxii, sdxii, + smeta, bgeta, aeta, veta, sdeta, ipage, page) call diagrm (lu, nconv, xsmall, xbig, ysmall, ybig) if (nconv .eq. 0) then do 1007 i = 1, 2 write (lu,*) ' ALL of your locations are out of', + ' bounds.' lu = nout 1007 continue endif elseif (itype .eq. 1) then * for file format types 1 and 2 lu = nout call report (lu, smxii, bgxii, axii, vxii, sdxii, + smeta, bgeta, aeta, veta, sdeta, ipage, page) call diagrm (lu, nconv, xsmall, xbig, ysmall, ybig) elseif (itype .eq. 2) then * itype = 2, (free format input, free format output) does not have * anything written to the output file in this subroutine elseif (itype .eq. 3 .or. itype .eq. 4 .or. + itype .eq. 5) then * itype = 3, itype = 4, and itype = 5 the ngs horizontal blue book * file formats do not have a * report written to the output file endif endif return end subroutine fgrid (xpt, ypt, dx, dy, xmax, xmin, + ymax, ymin, xgrid, ygrid, irow, jcol, nogo) ********************************************************************** ** subroutine fgrid: identifies the local grid square for interp. * ********************************************************************** * this subroutine is designed to identify the grid square in which a * particular point is located and get the corner coordinates * converted into the index coordinate system. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) double precision xpt, ypt real dx, dy, xmax, xmin, ymax, ymin, xgrid, ygrid integer irow, jcol logical nogo nogo = .false. * check to see it the point is outside the area of the gridded data if (xpt .ge. dble(xmax) .or. xpt .le. dble(xmin) .or. + ypt .ge. dble(ymax) .or. ypt .le. dble(ymin) ) then nogo = .true. * write (*,*) '***The point is out of bounds***' goto 200 endif * calculate the coordinate values for the point to be interpolated * in terms of grid indices xgrid = sngl(xpt - dble(xmin) )/dx + 1.e0 ygrid = sngl(ypt - dble(ymin) )/dy + 1.e0 * find the i,j values for the sw corner of the local square irow = ifix(ygrid) jcol = ifix(xgrid) 200 return end subroutine fhelp *** print information about the formats of the input data *** file types used by deflec. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea, itype, ios parameter (mxarea =18) character*1 ans integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) ***************************** * chose the input file format ***************************** 9001 write (luout,*) ' ' write (luout,*) ' What format do you want information about?' write (luout,*) ' 1) Free Format Type 1' write (luout,*) ' 2) Free Format Type 2' write (luout,*) ' 3) NGS Blue Book (old) format (default)' write (luout,*) ' 4) NGS Blue Book (new) format' write (luout,*) ' 5) NGS Blue Book (extended) format' read (luin,'(a1)') ans if (ans .eq. ' ') then itype = 3 else read (ans,90,err=9940,iostat=ios) itype 90 format (i1) if (itype .gt. 5 .or. itype .lt. 1) then * not a good answer - try again. write (luout,*) ' Gotta pick one of these -', + ' sorry try again.' goto 9001 endif endif * print information if (itype .eq. 1) then ******************************* * for file format itype = 1 * free format type1 input file * pretty output format ******************************* write (luout,2) * 2 format ('1') 2 format (' ') write (luout, 110) 110 format (' Free Format Type 1 - The first 40 characters of', + ' the input data record may', /, + ' contain the station name or be blank. The rest of', + ' the record (columns 41-80)', /, + ' must contain the latitude and longitude. They may', + ' be given in (1) decimal', /, + ' degrees; (2) integer degrees and decimal minutes,', + ' or (3) integer degrees,', /, + ' integer minutes, and decimal seconds. the decimal', + ' portion of the latitude', /, + ' MUST contain a decimal point as it is used to', + ' determine which is the last', /, + ' number forming part of the latitude. The output', + ' will be in "pretty" format.', /) write (luout, 120) 120 format (' The following three records are examples of valid', + ' input records:', /, + ' <------------ Columns 1-40 ------------>', + '<------------ Columns 41-80----------->', /, + ' AAA 34.', + '4444444 98.8888888', /, + ' BBB 25', + ' 55.55555 76 56.66666', /, + ' CCC 45 45', + ' 45.555 111 11 11.111', /) write (luout,931) 931 format (15x, ' (Hit RETURN to continue.)') read (luin,'(a1)') ans write (luout,2) write (luout, 130) 130 format (' The following is an example of the output. Note', + ' that with Free Format Type 1', /, + ' data, both the input latitude and the', + ' longitude are expressed in', /, + ' degrees, minutes, and seconds regardless of the', + ' method of input.', /) write (luout, 140) 140 format (' Computation #: 1', + ' Region: deflec96SCU(NGS)', //, + ' Station Name: AAA', //, + ' Latitude ', + 'Longitude Xi Deflec.', + ' Eta Deflec. Hor. Laplace', /, + ' ', + ' (arc sec) (arc sec) (arc sec)', //, + ' 34 26 39.99984 98 53', + ' 19.99968 -10.90 -5.50 3.77', /) elseif (itype .eq. 2) then ******************************* * for file format itype = 2 * free format type2 input file * free format output ******************************* write (luout,2) write (luout, 210) 210 format (' Free Format Type 2 - The first 32 characters of', + ' the input data record must', /, + ' contain the latitude and longitude. They may be', + ' given in (1) decimal degrees;', /, + ' (2) integer degrees and decimal minutes, or (3)', + ' integer degrees, integer', /, + ' minutes, and decimal seconds. The decimal portion', + ' of the latitude must contain', /, + ' a decimal point as it is used to determine which', + ' is the last number forming', /, + ' part of the latitude. The next 24 characters', + ' (columns 33-56) will hold the', /, + ' calculated deflections (Xi, Eta) and the horizontal', + ' Laplace correction. The ', /, + ' rest of the input record (columns 57-80) may contain', + ' the station name or be ', /, + ' blank. The output will be in the same format as the', + ' input but will contain ', /, + ' the deflections and Laplace corrections.', /) write (luout,931) read (luin,'(a1)') ans write (luout,2) write (luout, 220) 220 format (' The following three records are examples of', + ' valid input records:', //, + ' <------------ Columns 1-40 ------------>', + '<------------ Columns 41-80----------->', /, + ' 45 45 45.55555 111 11 11.11111 ', + 'one', /, + ' 25 55.5555555 76 56.6666666 ', + 'two', /, + ' 34.444444444 98.888888888 ', + 'three') write (luout, 230) 230 format (/, ' The following is an example of the output:', //, + ' DEFLEC Version 3.00 - NGS deflec96 Model', /, + ' <-------- Columns 1-32 -------->', + '<---- Columns 33-56---->', + '<---- Columns 57-80--->', /, + ' 45 45 45.55555 111 11 11.11111 6.08 -1.40', + ' 1.44one', /, + ' 25 55.5555555 76 56.6666666 2.98 22.22', + ' -10.80two', /, + ' 34.444444444 98.888888888 -10.90 -5.50', + ' 3.77three', /) elseif (itype .eq. 3) then **************************************** * for input file itype = 3 * the horizontal blue book specification **************************************** * the following two format statements logically belong together * but are split in order to have fewer that 19 continuation lines write (luout,2) write (luout, 310) 310 format (' NGS Horizontal ''Old'' Blue Book format - *80*', + ' (Control Point) and *85* (Defle', /, + 'ction) records. DEFLEC 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 DEFLEC - the *80*', + ' records and all other records are') write (luout, 311) 311 format (' passed through without change to the output file.', + ' On the *85* records, the SSN', /, + ' is used to determine which is the corresponding', + ' the *80* record. The', /, + ' information to the right of the SSN (columns', + ' 41-80) in the *85* record is', /, + ' overwritten by DEFLEC. On the *80* records, the', + ' direction of the latitude must', /, + ' be north positive (''N'') and the direction of the', + ' longitude must be west', /, + ' positive (''W'').', /) write (luout, 320) 320 format (' This format is out of date and should no longer be', + ' used for submitting data', /, + ' to the National Geodetic Survey.', //, + ' This format was defined in:', /, + ' ''Input Formats and Specifications of the', + ' National Geodetic Survey Data Base''', /, + ' ''Volume 1. Horizontal Control Data''.', /, + ' Published by the Federal Geodetic Control Committee', + ' in 1980 and available', /, + ' from: the National Geodetic Survey, NOAA,', + 'Silver Spring, MD 20910.', /) write (luout,931) read (luin,'(a1)') ans write (luout,2) write (luout, 330) 330 format (' The following input example is *80* and *85* records', + ' from an ''old'' Blue Book', /, + ' file:', //, + '004560*80*096 KNOXVILLE CT HSE', + ' 411906578 N0930548534 W 277 MIA33', /, + '004560*85*096 ABC bogus deflection ', + ' 10S999 280E999') write (luout, 340) 340 format (/, ' The following example is of the output *85*', + ' record with the new, interpolated', /, + ' deflections.', //, + '004560*85*096 NGS program DEFLEC deflc96NCL(NGS)', + ' 15N100 325E100',/) elseif (itype .eq. 4) then * the following two format statements logically belong together * but are split in order to have fewer that 19 continuation lines write (luout,2) write (luout,410) 410 format (' NGS Horizontal ''New'' Blue Book format - *80*', + ' (Control Point) and *85* (Defle', /, + 'ction) records. DEFLEC 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 DEFLEC - the *80*', + ' records and all other records are') write (luout, 411) 411 format (' passed through without change to the output file.', + ' on the *85* records, the SSN', /, + ' is used to determine which is the corresponding', + ' the *80* record. The', /, + ' information to the right of the SSN (columns', + ' 41-80) in the *85* record is', /, + ' overwritten by DEFLEC. On the *80* records, the', + ' direction of the latitude must', /, + ' be north positive (''N'') and the direction of the', + ' longitude must be west', /, + ' positive (''W'').', /) write (luout, 420) 420 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 January, 1989 and', /, + ' available from: the National Geodetic Survey, NOAA,', + ' Silver Spring, MD 20910.', /) write (luout,931) read (luin,'(a1)') ans write (luout,2) write (luout, 430) 430 format (' The following input example is *80* and *85* records', + ' from a ''new'' Blue Book', /, + ' file:', //, + '004560*80*0096KNOXVILLE CT HSE ', + '411906578 N0930548534 W 277 MIA33', /, + '004560*85*0096ABC bogus deflection ', + ' 10S999 280E999') write (luout, 440) 440 format (/, ' The following example is of the output *85*', + ' record with the new, interpolated', /, + ' deflection.', //, + '004560*85*0096NGS program DEFLEC deflc96NCL(NGS)', + ' 15N100 325E100',/) elseif (itype .eq. 5) then write (luout,2) write (luout,510) 510 format (' NGS Horizontal ''Extended'' Blue Book format', + ' (BIGADJUST format) - *80* (Control', /, + ' Point) and *85* (Deflection) records. Only the', + ' *80* records are used from', /, + ' the input file. The output consists only of *85*' + ' records. The SSN''s for the', /, + ' *85* records are taken from the *80* records - that', + ' is, there is an *85*', /, + ' record for each *80* record, regardless of any', + ' *85* records in the input file.', //, + ' This format is for internal NGS use only!!!', /) endif write (luout,*) ' Do you want more information (Y/N)?' write (luout,*) ' (Default is Y)' read (luin,'(a1)') ans if (ans .ne. 'n' .and. ans .ne. 'N') goto 9001 return * error message 9940 write (luout,*) ' Gotta pick ''1'' or ''2'' or ''3'' -', + ' sorry try again.' goto 9001 end subroutine getpt (nconv, itype, name, idla, imla, sla, idlo, + imlo, slo, xpt, ypt, eof, nopt) * get the name, latitude, and longitude of a point either interactively * or from an input data file * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) double precision xpt, ypt real sla, slo integer nconv, itype integer idla, imla, idlo, imlo character*80 name character*1 ans logical eof, nopt integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) eof = .false. nopt = .false. if (itype .eq. 0) then ************************************* * for interactive use - no input file ************************************* if (nconv .ge. 1) then write (luout,*) ' Do you want to do another deflection', + ' computation (Y/N)?' write (luout,*) ' (Default is Y)' read (luin,'(a1)') ans if (ans .eq. 'n' .or. ans .eq. 'N') goto 9999 endif * get a point (x,y) to compute call askpt (nconv, name, idla, imla, sla, + idlo, imlo, slo, xpt, ypt, eof, nopt) if (nopt) goto 9000 elseif (itype .eq. 1) then * free format type 1 call type1 (name, idla, imla, sla, idlo, imlo, slo, + xpt, ypt, eof, nopt) elseif (itype .eq. 2) then * free format type 2 call type2 (idla, imla, sla, idlo, imlo, slo, + xpt, ypt, eof, nopt) elseif (itype .eq. 3 .or. itype .eq. 4 .or. + itype .eq. 5 ) then * itype=3 is ngs old horizontal blue book * itype=3 is ngs new horizontal blue book * itype=5 is ngs extended horizontal blue book call type3 (itype, name, idla, imla, sla, idlo, imlo, slo, + xpt, ypt, eof, nopt) endif ********************************************************* * change the longitude to positive east for interpolation ********************************************************* xpt = -xpt 9000 return * end of file 9999 continue eof = .true. goto 9000 end subroutine grids * this subroutine opens the deflec grids that are requested in * a file named 'area.par' (if it exists). * area.par (or AREA.PAR) will be read for the names * and locations of the gridded deflection data set files. * the data in the area.par file has the following format: * columns 1-15 contain the name of the area. this name may * contain blanks or any other characters. * columns 16-80 (the rest of the record) contain the base name of the * grid files. that is the '.xii' and '.eta' extensions are not * included - it is added by this subroutine before each file is * opened. * comments are indicated by an '*' in column 1, blank records are * also treated as comments. comment records are printed to the * output file but otherwise ignored. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) character*80 b80 character*20 b20 parameter (b20 = ' ', b80 = b20//b20//b20//b20) real dx1, dy1, xmax1, xmin1, ymax1, ymin1 integer nr1, nc1 integer ios, i, iflag1, iflag2, n1, n2, leng, ierr, itemp character*80 card, ccard, dum character*65 afile, gfile character*15 aarea logical nogo, gflag character*15 areas common /areas/ areas(mxarea) integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea data dum / b80 / data iflag1 /1/, iflag2 /2/ *********************************************************************** * try and open the 'area.par' file, if not then try and open 'area.par' *********************************************************************** gfile = 'area.par' open (napar,file=gfile,form='formatted',status='old', + access='sequential',err=9100,iostat=ios) * file containing grid names exits 10 gflag = .true. write (luout, 80) 80 format ( '---> Data grids named in the AREA.PAR file <---', /, + ' # AREA NAME LOCATION', /, 1x, 79('=') ) do 120 i = 1, mxarea 100 read (napar,110,err=9000,end=9000,iostat=ios) card 110 format (a80) * check for comment records and blank records if ( card(1:1) .eq. '*' ) then call nblank (card, iflag2, n2) write (luout,'(5x, a)') card(1:n2) goto 100 elseif ( card .eq. b80 ) then ***** write (luout,*) ' ' goto 100 endif * get area name and basename of file (i.e. location without extensions) dum = card aarea = dum(1:15) call nblank (card(16:80), iflag1, n1) dum(1:65) = card(15+n1:80) leng = 65 afile = ccard(dum, leng, ierr) if (ierr .ne. 0) STOP 'Coding error in grids -- AFILE' if (afile .eq. ' ') goto 9000 * try and open files itemp = narea + 1 call openfl (afile, itemp, nogo, dx1, dy1, + xmax1, xmin1, ymax1, ymin1, nr1, nc1, card, gflag) if (.not. nogo) then * files opened ok and variables read narea = itemp areas(narea) = aarea *** patch for inexact headers (due to 2' spacing) 8-oct-96 dgm idx1=idnint(dx1*3600.d0) dx(NAREA) = dble(idx1)/3600.d0 idy1=idnint(dy1*3600.d0) dy(NAREA) = dble(idy1)/3600.d0 ***** dx(narea) = dx1 ***** old code ***** dy(narea) = dy1 ***** old code xmax(narea) = xmax1 xmin(narea) = xmin1 ymax(narea) = ymax1 ymin(narea) = ymin1 nr(narea) = nr1 nc(narea) = nc1 call nblank (card, iflag2, n2) write (luout,140) narea, card(1:n2) 140 format (2x, i2, 2x, a) endif 120 continue 9000 return * 'area.par' does not exist, try the name 'AREA.PAR' * if it exists, open it and continue, if it does not exist, return. 9100 continue gfile = 'AREA.PAR' open (napar,file=gfile,form='formatted',status='old', + access='sequential',err=9000,iostat=ios) goto 10 end subroutine headr (vrsion) *** this subroutine prints the header information and the disclaimer * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) real vrsion character*1 ans integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) write (luout,920) 920 format (10x, ' Welcome', /, + 10x, ' to the', /, + 10x, ' National Geodetic Survey', /, + 10x, ' Deflection Interpolation Program', //, + 10x, ' NGS DEFLEC96 Model', //, + 10x, ' (can also interpolate DMEX97,',/, + 10x, ' and DCAR97 models)',/) write (luout,930) vrsion 930 format ( /,10x, ' (Version', f5.2, ')', /, + 10x, ' July 16, 1997', /, + 10x, ' Dru Smith, Ph.D.',/, + 10x, ' Dennis Milbert, Ph.D.',/, + 10x, ' Systems Development Divison',/) write (luout,931) 931 format (10x, ' (Hit RETURN to continue.)') read (luin,'(a1)') ans write (luout,2) * 2 format ('1') 2 format (' ') write (luout,932) 932 format (/, ' 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 (luout,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.', /) write (luout,931) read (luin,'(a1)') ans write (luout,2) return end subroutine hms (dd, id, im, s) * use this to change from decimal degrees (double precision) * to integer degrees, integer minutes, and decimal seconds (real) * seconds are assumed to have no more than 5 decimal places * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) real small parameter (small = 1.e-5) double precision dd, temp real s integer id, im id = idint(dd) temp = ( dd - dble(id) )*60.0d0 im = idint(temp) s = sngl( ( temp - dble(im) )*60.0d0 ) if (im .eq. 60) then im = 0 id = id + 1 endif if (s .lt. small) s = 0.e0 if (s .gt. (60.e0-small) ) then s = 0.e0 im = im + 1 endif return end subroutine initl (screen, page, name, ipage, itype, + smxii, bgxii, axii, vxii, sdxii, sdxii2, + smeta, bgeta, aeta, veta, sdeta, sdeta2, + xsmall, xbig, ysmall, ybig) *** this subroutine initializes all the variables needed in deflec * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) character*20 b20 character*80 b80 parameter (b20 = ' ', b80 = b20//b20//b20//b20) double precision xsmall, xbig, ysmall, ybig double precision axii, vxii double precision sdxii, sdxii2 real smxii, bgxii double precision aeta, veta double precision sdeta, sdeta2 real smeta, bgeta integer ipage, itype character*80 name logical page, screen character*80 card common /curnt/ card integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) * initialize card variable in common curnt to blank card = b80 * set the logical units for input/output common inout * note that the unit numbers for the data grids (luarea) are defined * in subroutine openfl as the numbers 11 through 2*mxarea luin = 5 luout = 6 nout = 101 nin = 102 napar = 103 ****************************************************************** * initialize ****************************************************************** * defaults: screen = .true. => send results to screen * page = .false. => don't start a new page in the output file * name = ' ' => no station name * ipage = 0 => current output file page number is 0 * itype = 0 => interactive input of points screen = .true. page = .false. name = ' ' ipage = 0 itype = 0 smxii = 1.0e10 bgxii = -1.0e10 axii = 0.0d0 vxii = 0.0d0 sdxii = 0.0d0 sdxii2 = 0.0d0 smeta = 1.0e10 bgeta = -1.0e10 aeta = 0.0d0 veta = 0.0d0 sdeta = 0.0d0 sdeta2 = 0.0d0 xsmall = 1.0d10 xbig = -1.0d10 ysmall = 1.0d10 ybig = -1.0d10 return end subroutine interp (nogo, resp, xpt, ypt, xii, eta, itype) *** compute biquadratic interpolation, floating point version *** logic not tested for '0/360 meridian wraparound' * this code is adapted from that of dennis milbert - 6/5/90 * dennis(xpt,ypt,dx,dy,xmax,xmin,ymax,ymin, * & zee,nr,nc,luarea) * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mncbin, mxarea parameter (mncbin = 1000, mxarea =18) character*20 b20 character*80 b80 parameter (b20 = ' ', b80 = b20//b20//b20//b20) double precision xpt, ypt real buf, xii, eta real x, y, fx0, fx1, fx2, qterp2 real dx0, dy0, xmax0, ymax0, xmin0, ymin0 real xgrid, ygrid, val real tee1, tee2, tee3, tee4, tee5, tee6, tee7, tee8, tee9 integer nr0, nc0 integer iarea, irow, jcol, ifile1, ifile2, idum, i, j, idim integer ix, ix1, ix2, jy, jy1, jy2 integer iflag1, iflag2, n1, n2, itype character*15 resp logical nogo, flag dimension buf(mncbin) character*15 areas common /areas/ areas(mxarea) character*80 card common /curnt/ card integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea data iflag1 /1/, iflag2 /2/, flag /.false./ ****************************************************************** * initialize ****************************************************************** nogo = .false. **************************************************** * read where to get the data and how it is organized **************************************************** * check to see which set of gridded files xpt,ypt is in. do 100 iarea = 1, narea dx0 = dx(iarea) dy0 = dy(iarea) xmax0 = xmax(iarea) xmin0 = xmin(iarea) ymax0 = ymax(iarea) ymin0 = ymin(iarea) nr0 = nr(iarea) nc0 = nc(iarea) call fgrid (xpt, ypt, dx0, dy0, xmax0, xmin0, + ymax0, ymin0, xgrid, ygrid, irow, jcol, nogo) if (.not. nogo) goto 200 100 continue * not in any of the grid areas nogo = .true. goto 950 200 continue * point in area number iarea and named areas(iarea) ifile1= luarea(2*iarea-1) ifile2= luarea(2*iarea) resp = areas(iarea) if (ypt .lt. ymin0 .or. ypt .gt. ymax0 .or. + xpt .lt. xmin0 .or. xpt .gt. xmax0) then val=1.d30 else *** within grid boundaries, get indices in the grid ix = (xpt - xmin0)/dx0 + 1.d0 ix1 = ix + 1 ix2 = ix + 2 *** check if edge collision 1 if (ix2 .gt. nc0) then ix = ix - 1 ix1 = ix1 - 1 ix2 = ix2 - 1 go to 1 endif x = (xpt - dx0*(ix - 1) - xmin0)/dx0 *** move grid to get nearer center of 3 x 3 if (x .lt. 0.5 .and. ix .gt. 1) then ix = ix - 1 ix1 = ix1 - 1 ix2 = ix2 - 1 x = x + 1. endif if ( x .lt. 0. .or. x .gt. 2.) stop 'interp - 1' *** now do it for y jy = (ypt - ymin0)/dy0 + 1.d0 jy1 = jy + 1 jy2 = jy + 2 2 if (jy2 .gt. nr0) then jy = jy - 1 jy1 = jy1 - 1 jy2 = jy2 - 1 go to 2 endif y = (ypt - dy0*(jy - 1) - ymin0)/dy0 if (y .lt. 0.5 .and. jy .gt. 1) then jy = jy - 1 jy1 = jy1 - 1 jy2 = jy2 - 1 y = y + 1. endif if (y .lt. 0. .or. y .gt. 2.) stop 'interp - 2' *** get the nine data points needed from gridded deflection files !!! **************** * xii deflection **************** *** lower boundary read (ifile1,rec=jy+1) idum, (buf(j), j=1,nc0) tee1 = buf(ix) tee2 = buf(ix+1) tee3 = buf(ix+2) *** middle row read (ifile1,rec=jy1+1) idum, (buf(j), j=1,nc0) tee4 = buf(ix) tee5 = buf(ix+1) tee6 = buf(ix+2) *** upper boundary read (ifile1,rec=jy2+1) idim, (buf(j), j=1,nc0) tee7 = buf(ix) tee8 = buf(ix+1) tee9 = buf(ix+2) *** compute the value fx0 = qterp2(x, tee1, tee2, tee3) fx1 = qterp2(x, tee4, tee5, tee6) fx2 = qterp2(x, tee7, tee8, tee9) xii = qterp2(y, fx0, fx1, fx2) **************** * eta deflection **************** *** lower boundary read (ifile2,rec=jy+1) idum, (buf(j), j=1,nc0) tee1 = buf(ix) tee2 = buf(ix+1) tee3 = buf(ix+2) *** middle row read (ifile2,rec=jy1+1) idum, (buf(j), j=1,nc0) tee4 = buf(ix) tee5 = buf(ix+1) tee6 = buf(ix+2) *** upper boundary read (ifile2,rec=jy2+1) idim, (buf(j), j=1,nc0) tee7 = buf(ix) tee8 = buf(ix+1) tee9 = buf(ix+2) *** compute the value fx0 = qterp2(x, tee1, tee2, tee3) fx1 = qterp2(x, tee4, tee5, tee6) fx2 = qterp2(x, tee7, tee8, tee9) eta = qterp2(y, fx0, fx1, fx2) endif 9999 return * error messages 950 continue if (itype .ne. 0) then call nblank (card, iflag1, n1) call nblank (card, iflag2, n2) write (luout,955) card(n1:n2) 955 format (' *** This point is out of bounds ***', /, + 1x, '''', a, '''') else write (luout,960) 960 format (' *** The point is out of bounds ***') endif * write out grid areas for the first out-of-bounds error message if (.not.flag .or. itype .eq. 0) then write (luout,*) ' It must be within one of the following grid', + ' areas;' write (luout,975) 975 format (18x, 7x, 'Latitude', 7x, 'Longitude', /, + 5x, 'Area Name', 5x, 2(5x, 'Min', 4x, 'Max', 1x), + '(degrees)' ) do 970 i = 1, narea write (luout,965) areas(i), nint( ymin(i) ), nint( ymax(i) ), + nint( -xmax(i) ), nint( -xmin(i) ) 965 format (1x, '''', a15, '''', 2(2x, 2i7) ) 970 continue flag = .true. endif write (luout,*) ' ' goto 9999 end subroutine iparms (itype, screen) *** this subroutine interactively requests for information *** needed by deflec. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) character*80 b80 character*20 b20 parameter (b20 = ' ', b80 = b20//b20//b20//b20) character*80 infile, ofile character*1 ans integer itype integer iflag1, iflag2, n1, n2, ios logical screen, eflag integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) data iflag1 /1/, iflag2 /2/ ********************************** * get the name for the output file ********************************** 14 write (luout,*) ' What is the name of the file which will', + ' contain the program results?' write (luout,*) ' The default name is ''deflec.out''.' read (luin,'(a80)') ofile if (ofile .eq. b80) ofile = 'deflec.out' inquire (file=ofile, exist=eflag) if (eflag) then call nblank (ofile, iflag1, n1) call nblank (ofile, iflag2, n2) write (luout,*) ' The file ''', ofile(n1:n2), '''' write (luout,*) ' already exists. Do you want to overwrite', + ' it (Y/N)?' write (luout,*) ' (Default is Y)' read (luin, '(a1)') ans if (ans .eq. 'n' .or. ans .eq. 'N') goto 14 endif open (nout,file=ofile,form='formatted',status='unknown', + access='sequential',err=9910,iostat=ios) ********************************** * get the name for the input file ********************************** write (luout,*) ' ' 13 write (luout,*) ' Do you have an input data file (Y/N)?' write (luout,*) ' (Default is N)' read (luin,'(a1)') ans if (ans .eq. 'y' .or. ans .eq. 'Y') then write (luout,*) ' ' write (luout,*) ' What is the name of the input data file?' write (luout,*) ' The default name is ''BBOOK''.' read (luin,'(a80)') infile if (infile .eq. b80) infile = 'BBOOK' open (nin,file=infile,form='formatted',status='old', + access='sequential',err=9920,iostat=ios) ***************************** * chose the input file format ***************************** 9001 continue call nblank (infile, iflag2, n2) write (luout,*) ' ' write (luout,*) ' What is your file format?' write (luout,*) ' 0) Help - file format information' write (luout,*) ' 1) Free Format Type 1' write (luout,*) ' 2) Free Format Type 2' write (luout,*) ' 3) NGS Blue Book (old) format (default)' write (luout,*) ' 4) NGS Blue Book (new) format' write (luout,*) ' 5) NGS Blue Book (extended) format' read (luin,'(a1)') ans if (ans .eq. ' ') then * ngs horizontal blue book format itype = 3 else read (ans,347,err=9940,iostat=ios) itype 347 format (i1) if (itype .gt. 5 .or. itype .lt. 0) then * not a good answer - try again. write (luout,*) ' Gotta pick one of these -', + ' sorry try again.' goto 9001 endif endif * get help information if (itype .eq. 0) then call fhelp goto 9001 endif ******************************** * check for a screen output also ******************************** write (luout,*) ' ' write (luout,*) ' Do you want the results written to the', + ' terminal screen as well as to' write (luout,*) ' the output file (Y/N)?' write (luout,*) ' (Default is N)' read (luin,'(a1)') ans if (ans .ne. 'y' .and. ans .ne. 'Y') screen = .false. goto 9002 * error message 9940 write (luout,*) ' Gotta pick ''1'' or ''2'' or ''3'' -', + ' sorry try again.' goto 9001 9002 endif return * error message 9910 continue call nblank (ofile, iflag1, n1) call nblank (ofile, iflag2, n2) write (luout,9915) ios, ofile(n1:n2) 9915 format (' ERROR (', i5, ') - The operating system could not', + ' open the file ', /, + ' ''', a, '''', /, + ' Try again.', /) goto 14 9920 continue call nblank (infile, iflag1, n1) call nblank (infile, iflag2, n2) write (luout,9915) ios, infile(n1:n2) goto 13 end subroutine mloop (nconv, ipage, itype, vrsion, + xii, sdxii, sdxii2, smxii, bgxii, + eta, sdeta, sdeta2, smeta, bgeta, + xsmall, xbig, ysmall, ybig, + page, screen, name) ********************************************************************** * this subroutine loops through the input data (either an input data * * file or interactively), calculates the deflection values, * * updates the minimum, maximum, and statistical summations, and then * * prints the results to the output file and/or the screen. * ********************************************************************** * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) double precision xsmall, xbig, ysmall, ybig, xpt, ypt real vrsion, sla, slo double precision sdxii, sdxii2 real xii, smxii, bgxii, xii2 double precision sdeta, sdeta2 real eta, smeta, bgeta, eta2 integer nconv, ipage, itype integer idla, imla, idlo, imlo character*80 name character*15 resp logical page, nogo, screen, nopt, eof real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) * set defaults for those variables not used by every format type ******************************************************************** * begin the computation loop for each interpolation * do until end of file or no more computations requested ******************************************************************** nconv = 0 160 continue page = .false. ******************************************** * get the name and location of another point ******************************************** call getpt (nconv, itype, name, idla, imla, sla, idlo, + imlo, slo, xpt, ypt, eof, nopt) if (nopt) goto 155 if (eof) goto 9999 ********************** * do the interpolation ********************** nogo = .false. call interp (nogo, resp, xpt, ypt, xii, eta, itype) **************************************************** * check to see if this point is valid * if nogo is true then get another point and don't * do the computation - point is out of bounds * if nogo is not true then proceed - estimate made **************************************************** if (nogo) goto 155 nconv = nconv + 1 ******************************************** * change the longitude back to positive west * and change back to d.m.s format ******************************************** xpt = -xpt ************************** * do the little statistics ************************** * first, the basics * arc sec.... xii2 = ( dble(xii) )**2 sdxii2 = sdxii2 + xii2 sdxii = sdxii + dble(xii) eta2 = ( dble(eta) )**2 sdeta2 = sdeta2 + eta2 sdeta = sdeta + dble(eta) * then the ranges xsmall = dmin1( xsmall, xpt) xbig = dmax1( xbig , xpt) ysmall = dmin1( ysmall, ypt) ybig = dmax1( ybig , ypt) smxii = min( smxii, xii) bgxii = max( bgxii, xii) smeta = min( smeta, eta) bgeta = max( bgeta, eta) ********************************** ** write to output file and screen ********************************** call wrtpt (itype, nconv, vrsion, name, + idla, imla, sla, idlo, imlo, slo, + xii, eta, resp, ipage, page, screen) ********************** * start the loop again ********************** 155 goto 160 9999 return end subroutine nblank (a, iflag, nblk) *** return position of last non-blank in string (iflag = 2) *** or position of first non-blank in string (iflag = 1) * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer iflag, nblk, leng, iblk character*(*) a leng = len(a) if (iflag .eq. 2) then do 1 iblk = leng, 1, -1 if ( a(iblk:iblk) .ne. ' ' ) then nblk = iblk return endif 1 continue elseif (iflag .eq. 1) then do 2 iblk = 1, leng, +1 if ( a(iblk:iblk) .ne. ' ' ) then nblk = iblk return endif 2 continue endif * string contains all blanks nblk = 0 return end subroutine ngrids (nodata) * this subroutine opens the deflec grids. * two file are necessary for each area. * if a file named area.par exists it will be read for the names and * locations of the gridded data. the format of the data in * the area.par file is given in the grids subroutine. * if the area.par file does not exist, or there is still room in the * arrays in the gdinfo common then the default area names used. * implicit real (a-h, o-z) * implicit integer (i-n) * implicit undefined (a-z) integer mxarea parameter (mxarea =18) character*1 ans logical nodata integer luin, luout, nout, nin, napar, luarea common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea) real dx, dy, xmax, xmin, ymax, ymin integer nr, nc, narea common /gdinfo/ dx(mxarea), dy(mxarea), xmax(mxarea), + xmin(mxarea), ymax(mxarea), ymin(mxarea), + nr(mxarea), nc(mxarea), narea *** dmex97 and dcar97 flags logical ldmex,ldcar common/scienc/ldmex,ldcar * initialize nodata = .false. narea = 0 write (luout,100) 100 format (' DEFLEC is now opening the files containing the', + ' gridded data.') * try to open the 'area.par' file in the subroutine grids call grids * if narea>=mxarea, then skip the section that opens the default files. * if narea