program defgrd *********************************************************************** * * * purposes: 1) to extract deflec subgrids from larger deflec grids. * * * * 2) to translate deflec grid files to binary from the * * ascii transfer format and visa versa. * * * * 3) to translate deflec grid files from binary or the * * ascii transfer format to the ascii graphics format. * * * * 4) to print information about deflec grids to a file. * * * * the program deflec reads binary grids of deflection * * of the vertical in the meridian and prime vertical. * * the first record in a grid file consists of header * * information. all the other records consist of real*4 * * grid numbers. the grid files are unformatted and * * direct access. * * * * deflec grids come in pairs and consist of the xi * * and eta deflection in seconds of arc. * * the names of the grids consist of * * the location of the grid (e.g. 'd90se' or 'd90nc') * * plus an extension. the extension for the grid of * * xi (meridian n/s) is '.xii' and the extension for the * * grid of eta (prime vertical e/w) is '.eta'. * * * * the ascii transfer files are used for portability. * * the first two records in the ascii file contain the * * header information. the rest of the records are * * formatted as 6f12.6. * * * * this program is intended to be used with deflec * * development. * * * * variables in this code that end in an 'a' are * * associated with a xi grid while variables that * * end in an 'o' are associated with an eta grid. * * * * all grids start and finish at whole degrees. * * note that there must be minimum number of degrees * * between the minimum and maximum longitude. this is * * explained further in subroutine nlon. * * * * version code: 1.00 * * * * version date: september 5, 1991 * * * * author: dennis g. milbert, ph.d. * * advanced geodetic science branch * * national geodetic survey, nos, noaa * * rockville, md 20852 * * * * original authors of program nadgrd: * * alice r. drew * * senior geodesist, horizontal network branch * * warren t. dewhurst, ph.d. * * lieutenant commander, noaa * * national geodetic survey, nos, noaa * * rockville, md 20852 * *********************************************************************** *********************************************************************** * * * 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) * nmax is the maximum size of the buffer. it is at least as large * as the largest number of columns in the input grid integer nmax parameter (nmax = 1024) real vrsion parameter (vrsion = 1.00e0) real xmin, xmax, ymin, ymax integer ncfrst, nclast, nrfrst, nrlast integer kin, kout character*56 rident character*32 fouta, fouto character*28 basein character*8 pgm logical flag integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto * initialize the input/output unit numbers luin = 5 luout = 6 nina = 11 nino = 12 nouta = 21 nouto = 22 * header call headr (vrsion) * open the input and output files flag = .false. call ifiles (kin, flag, basein, rident, pgm) if (flag) goto 9000 call ofiles (kin, kout, fouta, fouto, basein, rident, pgm) * get the variables for the new grid files. call nparms (kout, xmin, xmax, ymin, ymax, + nrfrst, nrlast, ncfrst, nclast) * get maximum and minimum z values in new grid area call zextrm (kin, kout, ncfrst, nclast, nrfrst, nrlast) *** copy from input files to output files if (kout .ne. 0) then call tonew (kin, kout, xmin, xmax, ymin, ymax, + ncfrst, nclast, nrfrst, nrlast, fouta, fouto) endif * close all files 9000 continue close (nina,status='keep') close (nino,status='keep') if (.not. flag) then close (nouta,status='keep') if (kout .ne. 0) close (nouto,status='keep') endif write (luout,9010) 9010 format (' End of program DEFGRD') stop end subroutine headr (vrsion) *** this subroutine prints the header information implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) real vrsion character*1 ans integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto write (luout,920) 920 format (12x, ' Welcome', /, + 12x, ' to the', /, + 12x, ' National Geodetic Survey', /, + 12x, ' Deflection Grid Manipulation Program.', + //, 12x, ' For use when DEFLEC Grids', /, + 12x, ' need to be translated', /, + 12x, ' between ASCII and binary', /, + 12x, ' or', /, + 12x, ' when grids covering a smaller areal extent', /, + 12x, ' need to be extracted from standard', /, + 12x, ' DEFLEC grids.', /, + 12x, ' or', /, + 12x, ' for obtaining information', /, + 12x, ' about the DEFLEC grids.') write (luout,930) vrsion 930 format (/, 12x, ' (version', f5.2, ')', /, + 12x, ' September 5, 1991', /, + 12x, ' Dennis G. Milbert, Ph.D.', /, + 12x, ' Advanced Geodetic Science Branch',/) write (luout,931) 931 format (12x, ' (Hit RETURN to continue.)') read (luin,'(a1)') ans write (luout,2) * 2 format ('1') 2 format (' ') write (luout,932) 932 format ( /, 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 (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 return end subroutine ifiles (kin, flag, basein, rident, pgm) * interactively get input file basenames. * open the input files implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) character*32 b32 parameter (b32 = ' ') real xmin0a, ymin0a, dxa, dya real angle, anglea integer nz, lrecl, kin integer iflag2, n2 integer nca, nra, nza character*56 rident character*32 fina, fino character*28 basein character*8 pgm character*1 ans logical flag integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nc, nr common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc data iflag2 /2/ * find out about the input files write (luout,2) * 2 format ('1') 2 format (' ') 220 write (luout,225) 225 format ( ' For the input files enter:', + /, ' ''A'' for ASCII transfer format.', + /, ' ''B'' for binary - DEFLEC data file format).', + /, ' (Default is B)') read (luin,'(a1)') ans if (ans .eq. ' ' .or. ans .eq. 'b' .or. ans .eq. 'B') then kin = 1 elseif (ans .eq. 'a' .or. ans .eq. 'A') then kin = -1 else write (luout, 210) ans 210 format (/, ' ERROR - ''', a1, ''' is not a legal answer.') goto 220 endif * get basename if (kin .eq. 1) then write (luout,110) 110 format ( ' Enter the basename for the pair of input DEFLEC', + ' grid files from which the', + /, ' new grids will be extracted. The ''.xii''', + ' and ''.eta'' extensions will be', + /, ' added to the basename by DEFGRD. The default', + ' basename is ''d90sw''.') elseif (kin .eq. -1) then write (luout,115) 115 format ( ' Enter the basename for the pair of input DEFLEC', + ' grid files from which the', + /, ' new grids will be extracted. The ''.xia''', + ' and ''.eaa'' extensions will be', + /, ' added to the basename by DEFGRD. The default', + ' basename is ''d90sw''.') endif read (luin,'(a28)') basein call nblank (basein, iflag2, n2) if (n2 .eq. 0) then basein = 'd90sw' n2 = 5 endif fina = b32 fina(1:n2) = basein(1:n2) fino = b32 fino(1:n2) = basein(1:n2) if (kin .eq. 1) then * open binary input files n2 = n2 + 4 fina(n2-3:n2) = '.xii' fino(n2-3:n2) = '.eta' lrecl = 256 open (nina,file=fina,form='unformatted',access='direct', + recl=lrecl,status='old',err= 900) open (nino,file=fino,form='unformatted',access='direct', + recl=lrecl,status='old',err= 900) read (nina,rec=1) rident, pgm, nca, nra, nza, xmin0a, dxa, + ymin0a, dya, anglea read (nino,rec=1) rident, pgm, nc, nr, nz, xmin0, dx, + ymin0, dy, angle elseif (kin .eq. -1) then * open ascii input files n2 = n2 + 4 fina(n2-3:n2) = '.xia' fino(n2-3:n2) = '.eaa' open (nina,file=fina,form='formatted',access='sequential', + status='old',err=900) open (nino,file=fino,form='formatted',access='sequential', + status='old',err=900) read (nina,140) rident, pgm read (nino,140) rident, pgm 140 format (a56, a8) read (nina, 150,err=900) nca, nra, nza, xmin0a, dxa, + ymin0a, dya, anglea read (nino, 150,err=900) nc, nr, nz, xmin0, dx, + ymin0, dy, angle 150 format (3i4, 5f12.5) endif * check for corrupted files if (nca .ne. nc .or. nra .ne. nr .or. nza .ne. nz .or. + xmin0a .ne. xmin0 .or. ymin0a .ne. ymin0 .or. + dxa .ne. dx .or. dya .ne. dy .or. anglea .ne. angle) + goto 960 if (nz .ne. 1 .or. angle .ne. 0.0) goto 960 * reopen binary files with correct record length if (kin .eq. 1) then close (nina) close (nino) lrecl = 4*(nc + 1) open (nina,file=fina,form='unformatted',access='direct', + recl=lrecl,status='old') open (nino,file=fino,form='unformatted',access='direct', + recl=lrecl,status='old') endif * write the (user-significant) input file variables to the screen xmax0 = dx * real( nc-1 ) + xmin0 ymax0 = dy * real( nr-1 ) + ymin0 if (kin .eq. 1) then write (luout,130) fina(1:n2), fino(1:n2) 130 format (/, ' Binary files ''', a, ''' and ''', a, ''' have', + ' been opened.') elseif (kin .eq. -1) then write (luout,135) fina(1:n2), fino(1:n2) 135 format (/, ' ASCII files ''', a, ''' and ''', a, ''' have', + ' been opened.') endif write (luout,180) ymin0, ymax0, -xmax0, -xmin0 180 format (/, ' Their minimum and maximum latitudes (+N) are: ', + 2f8.0, + /, ' Their minimum and maximum longitudes (+W) are: ', + 2f8.0, /) write (luout,200) 200 format (14x, '(Hit RETURN to continue.)') read (luin,'(a1)') ans write (luout,2) 999 return * error messages 900 write (luout,910) fina(1:n2), fino(1:n2) 910 format (/, ' *** ERROR *** ''', a, ''' and ''', a, ''' cannot', + ' be opened!', + /, ' Do you wish to try again (Y/N)?', + /, ' (Default is Y)') read (luin,'(a1)') ans if (ans .ne. 'n' .and. ans .ne. 'N') goto 220 flag = .true. goto 999 * corrupted input files 960 write (luout, 965) fina(1:n2), fino(1:n2) 965 format (/, 5x, ' *********** ERROR ***********', /, + ' The header information in grid files', /, + ' ''', a, ''' and ''', a, '''', /, + ' is incorrect! One or both of these files is', + ' corrupted or the wrong format.') flag = .true. goto 999 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 leng, iflag, iblk, nblk 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 nlat (kout, ymin, ymax) * get the new latitude variables implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * small is system dependent and should be a very small real number real small parameter (small = 1.e-5) real rcard real ymin, ymax integer leng, ierr, kout character*10 ans integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nr, nc common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc 110 write (luout,115) ymin0 115 format (/, ' Enter the minimum latitude, north positive.', + /, ' The default value is:', f8.0) read (luin,'(a10)') ans if (ans .eq. ' ') then ymin = ymin0 else leng = 10 ymin = rcard (ans, leng, ierr) endif write (luout,120) ymax0 120 format (/, ' Enter the maximum latitude, north positive.', + /, ' The default value is:', f8.0) read (luin,'(a10)') ans if (ans .eq. ' ') then ymax = ymax0 else leng = 10 ymax = rcard (ans, leng, ierr) endif * check values if (ymin .lt. ymin0 .or. ymin .gt. ymax0 .or. + ymax .lt. ymin0 .or. ymax .gt. ymax0) then write (luout, 130) ymin0, ymax0 130 format (/, ' *** ERROR *** ', + /, ' One or both of the values that you gave for the', + ' minimum and maximum', + /, ' latitudes are out of bounds! For the input', + ' DEFLEC grids', + /, ' the minimum and maximum latitudes (+N) are: ', + 2f8.0) goto 110 elseif ( ( ymin .gt. ymax ) .or. + ( kout .eq. -2 .and. ymin .eq. ymax) ) then write (luout, 140) ymin0, ymax0 140 format (/, ' *** ERROR *** ', + /, ' The minimum value that you gave for the latitude', + ' is less than the maximum!', + /, ' For the input DEFLEC grids', + /, ' the minimum and maximum latitudes (+N) are: ', + 2f8.0) goto 110 endif ********************************************************************* * the subgrid must be buffered by at least one whole degree of * the original grid - but must be wholy within the original grid. ********************************************************************* * modification note -- this is an affront to humanity -- DISABLED! ********************************************************************* * for ascii graphics format, buffer only to the nearest whole degree if (kout .ne. -2) then *** modified, dgm 5-sep-91 *** * ymin = aint(ymin - 1.e0 + small) * ymax = aint(ymax + 2.e0 - small) ymin = aint(ymin + small) ymax = aint(ymax + 1.e0 - small) else ymin = aint(ymin + small) ymax = aint(ymax + 1.e0 - small) endif if (ymin .lt. ymin0) ymin = ymin0 if (ymax .gt. ymax0) ymax = ymax0 return end subroutine nlon (kout, xmin, xmax) * get the new longitude variables. * the variables lmin is the minimum number of degrees between the * minimum longitude and the maximum longitude. the reason that * there must be this minimum is that in order to be read by deflec, * the extracted grid must have a logical record length of at least 96. * this is because the header record length is 96 bytes. * this variables is not used for the graphics output (kout=-2) * what happens is that the output file is opened with a record length * of 4*(nc2+1) and the first record written. the second record * is then written into the file starting at character 4*(nc2+1)+1. * where nc2 is the number of columns in the output file and is * determined by the longitude difference and the longitude increment. * since the values are written as real*4 numbers, there must be at * least 24 numbers in a row. * for example, the deflec files have dx=3'. defgrd requires that the * grid boundarys be whole degrees. in order to have at least 24 * number in a row, the (maximum-minimum) longitude must be at least 2 * degrees. implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * small is system dependent and should be a very small real number real small parameter (small = 1.e-5) real rcard real xmin, xmax integer leng, ierr, kout, lmin character*10 ans integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nr, nc common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc 160 write (luout,150) -xmax0 150 format (/, ' Enter the minimum longitude,', + ' west longitude positive.', + /, ' The default value is:', f8.0) read (luin,'(a10)') ans if (ans .eq. ' ') then xmax = xmax0 else leng = 10 xmax = rcard (ans, leng, ierr) xmax = -xmax endif if (kout .ne. -2) then lmin = int(24.e0*dx + small) if (lmin .lt. 1) lmin = 1 write (luout,155) lmin, -xmin0 155 format (/, ' Enter the maximum longitude,', + ' west longitude positive.', + /, ' The difference between the minimum and', + ' maximum longitudes', + /, ' must be at least', i3, ' degrees. The default', + ' value is:', f8.0) else write (luout,165) -xmin0 165 format (/, ' Enter the maximum longitude,', + ' west longitude positive.', + /, ' The default value is:', f8.0) endif read (luin,'(a10)') ans if (ans .eq. ' ') then xmin = xmin0 else leng = 10 xmin = rcard (ans, leng, ierr) xmin = -xmin endif * check values if (xmin .lt. xmin0 .or. xmin .gt. xmax0 .or. + xmax .lt. xmin0 .or. xmax .gt. xmax0) then write (luout, 170) -xmax0, -xmin0 170 format (/, ' *** ERROR *** ', + /, ' One or both of the values that you gave for the', + ' minimum and maximum', + /, ' longitudes are out of bounds! For the input', + ' DEFLEC grids', + /, ' the minimum and maximum longitudes (+W) are: ', + 2f8.0) goto 160 elseif ( ( xmin .gt. xmax ) .or. + ( kout .eq. -2 .and. xmin .eq. xmax) ) then write (luout, 180) -xmax0, -xmin0 180 format (/, ' *** ERROR *** ', + /, ' The minimum value that you gave for the longitude', + ' is less than the maximum!', + /, ' For the input DEFLEC grids', + /, ' the minimum and maximum longitudes (+W) are: ', + 2f8.0) goto 160 endif ********************************************************************* * the subgrid must be buffered by at least one whole degree of * the original grid - but must be wholy within the original grid. ********************************************************************* * modification note -- this is an affront to humanity -- DISABLED! ********************************************************************* * for ascii graphics format, buffer only to the nearest whole degree if (kout .ne. -2) then *** modified, dgm 5-sep-91 *** * xmin = aint(xmin - 2.e0 + small) * xmax = aint(xmax + 1.e0 - small) xmin = aint(xmin - 1.e0 + small) xmax = aint(xmax - small) else xmin = aint(xmin - 1.e0 + small) xmax = aint(xmax - small) endif if (xmin .lt. xmin0) xmin = xmin0 if (xmax .gt. xmax0) xmax = xmax0 * check range if ( kout .ne. -2 .and. (xmax-xmin) .lt. lmin ) then write (luout, 190) lmin, -xmax0, -xmin0 190 format (/, ' *** ERROR *** ', + /, ' The difference between the minimum and maximum', + ' longitudes must be at', + /, ' least', i3, ' degrees. ', + ' For the input DEFLEC grids', + /, ' the minimum and maximum longitudes (+W) are: ', + 2f8.0) goto 160 endif return end subroutine nparms (kout, xmin, xmax, ymin, ymax, + nrfrst, nrlast, ncfrst, nclast) * get the variables for the new grid files. implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * small is system dependent and should be a very small real number real small parameter (small = 1.e-5) real xmin, xmax, ymin, ymax integer kout, ncfrst, nclast, nrfrst, nrlast character*1 ans integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nc, nr common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc * if only printing grid information then set defaults and return if (kout .eq. 0) then ymin = ymin0 ymax = ymax0 xmax = xmax0 xmin = xmin0 nrfrst = 1 nrlast = nr ncfrst = 1 nclast = nc return endif *** get the area variables of the grid write (luout,2) * 2 format ('1') 2 format (' ') if (kout .ne. -2) then write (luout,100) 100 format ( ' The extracted DEFLEC grids will be automatically', + ' set to the', + /, ' whole degree by DEFGRD. Please enter', + ' the latitude and longitude', + /, ' extremes for your area of interest. Enter the', + ' values in decimal degrees.', + /, ' If you enter blanks, the extremes from the input', + ' files will be used.') else write (luout,105) 105 format ( ' The extracted DEFLEC grids will be automatically', + ' rounded to the', + /, ' whole degree by DEFGRD. Please enter the', + ' latitude and longitude', + /, ' extremes for your area of interest. Enter the', + ' values in decimal degrees.', + /, ' If you enter blanks, the extremes from the input', + ' files will be used.') endif * latitude 200 call nlat (kout, ymin, ymax) * longitude, change to east positive call nlon (kout, xmin, xmax) * ask if the new maximums and minimums are ok write (luout,2) if ( kout .ne. -2 .and. + ( xmin .ne. xmin0 .or. xmax .ne. xmax0 .or. + ymin .ne. ymin0 .or. ymax .ne. ymax0 ) ) then write (luout, 210) ymin, ymax, -xmax, -xmin 210 format ( ' For the output grids;', + /, ' the minimum and maximum latitudes ', + ' (+N) are: ', 2f8.0, + /, ' the minimum and maximum longitudes', + ' (+W) are: ', 2f8.0) else write (luout, 215) ymin, ymax, -xmax, -xmin 215 format ( ' For the output grids;', + /, ' the minimum and maximum latitudes ', + ' (+N) are: ', 2f8.0, + /, ' the minimum and maximum longitudes', + ' (+W) are: ', 2f8.0) endif write (luout,220) 220 format (/, ' If these values are OK, hit RETURN to continue.', + ' Enter any other character', + /, ' to re-enter the minimum and maximum latitudes', + ' and longitudes.') read (luin,'(a1)') ans if (ans .ne. ' ') goto 200 * calculate the new origin for the subgrid area, * calculate the number of columns and rows nrfrst = nint( (ymin - ymin0)/dy ) + 1 nrlast = nint( (ymax - ymin0)/dy ) + 1 ncfrst = nint( (xmin - xmin0)/dx ) + 1 nclast = nint( (xmax - xmin0)/dx ) + 1 return end subroutine ofiles (kin, kout, fouta, fouto, basein, rident, pgm) * interactively get output file basenames. implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) character*28 b28 parameter (b28 = ' ') character*32 b32 parameter (b32 = ' ') integer kin, kout integer iflag2, n2 character*56 rident character*32 fouta, fouto character*28 basein, basout, dbname character*8 pgm character*1 ans integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nc, nr common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc data iflag2 /2/ * find out about the output file(s) 220 write (luout,115) 115 format ( ' For the output files enter:', + /, ' ''A'' for ASCII transfer format.', + /, ' ''B'' for binary - DEFLEC data file format.', + /, ' ''G'' for ASCII graphics format. ', + ' There are five header information records,', + /, ' and the record lengths', + ' can be VERY large.' + /, ' ''I'' for input file information written to a', + ' single output file.', + /, ' (Default is B)') read (luin,'(a1)') ans if (ans .eq. ' ' .or. ans .eq. 'b' .or. ans .eq. 'B') then kout = 1 elseif (ans .eq. 'a' .or. ans .eq. 'A') then kout = -1 elseif (ans .eq. 'g' .or. ans .eq. 'G') then kout = -2 elseif (ans .eq. 'i' .or. ans .eq. 'I') then kout = 0 else write (luout, 210) ans 210 format (/, ' ERROR - ''', a1, ''' is not a legal answer.') goto 220 endif * default basename is the input basename unless the same format if (kout .ne. kin) then dbname = basein call nblank (dbname, iflag2, n2) else dbname = 'defgrd' n2 = 6 endif 180 if (kout .eq. 1) then * the extensions for the binary output files are .xii and .eta write (luout,160) dbname(1:n2) 160 format (/, ' Enter the file name for the new pair of', + ' DEFLEC grids. Enter only the', + /, ' basename. The ''.xii'' and ''.eta'' extensions', + ' will be added to the', + /, ' basename by DEFGRD. The default base file name', + ' is ''', a, '''.') read (luin,'(a28)') basout if (basout .eq. b28) basout = dbname if (kin .eq. 1 .and. basout .eq. basein) then write (luout,165) basout 165 format (/, ' ***** ERROR ****', + /, ' File ''', a, ''' is the name of the input file!', + /, ' You must choose another name.') goto 180 endif call nblank (basout, iflag2, n2) fouta = b32 fouta(1:n2) = basout(1:n2) fouto = b32 fouto(1:n2) = basout(1:n2) fouta(n2+1:n2+4) = '.xii' fouto(n2+1:n2+4) = '.eta' elseif (kout .eq. -1) then * the extensions for the ascii output files are .xia and .eaa write (luout,190) dbname(1:n2) 190 format (/, ' Enter the file name for the new pair of', + ' DEFLEC grids. Enter only the', + /, ' basename. The ''.xia'' and ''.eaa'' extensions', + ' will be added to the', + /, ' basename by DEFGRD. The default base file name', + ' is ''', a, '''.') read (luin,'(a28)') basout if (basout .eq. b28) basout = dbname if (kin .eq. -1 .and. basout .eq. basein) then write (luout,185) basout 185 format (/, ' ***** ERROR ****', + /, ' File ''', a, ''' is the name of the input file!', + /, ' You must choose another name.') goto 180 endif call nblank (basout, iflag2, n2) fouta = b32 fouta(1:n2) = basout(1:n2) fouto = b32 fouto(1:n2) = basout(1:n2) fouta(n2+1:n2+4) = '.xia' fouto(n2+1:n2+4) = '.eaa' elseif (kout .eq. -2) then * the extensions for the graphics output files are .lag and .log write (luout,170) dbname(1:n2) 170 format (/, ' Enter the file name for the new pair of', + ' DEFLEC grids. Enter only the', + /, ' basename. The ''.lag'' and ''.log'' extensions', + ' will be added to the', + /, ' basename by DEFGRD. The default base file name', + ' is ''', a, '''.') read (luin,'(a28)') basout if (basout .eq. b28) basout = dbname call nblank (basout, iflag2, n2) fouta = b32 fouta(1:n2) = basout(1:n2) fouta(n2+1:n2+4) = '.lag' fouto = b32 fouto(1:n2) = basout(1:n2) fouto(n2+1:n2+4) = '.log' elseif (kout .eq. 0) then * the extension for the information file is .inf write (luout,195) dbname(1:n2) 195 format (/, ' Enter the file name for the grid information', + ' file. Enter only the basename.', + /, ' The ''.inf'' extension will be added to the', + ' basename by DEFGRD. The default', + /, ' base file name is ''', a, '''.') read (luin,'(a28)') basout if (basout .eq. b28) basout = dbname call nblank (basout, iflag2, n2) fouta = b32 fouta(1:n2) = basout(1:n2) fouta(n2+1:n2+4) = '.inf' * write input header information to the information file open (nouta,file=fouta,form='formatted',access='sequential', + status='unknown') if (kin .eq. 1) then call nblank (basein, iflag2, n2) write (nouta,230) basein(1:n2)//'.xii', basein(1:n2)//'.eta' 230 format (/, ' Files ''', a, ''' and ''', a, + ''' are binary.') elseif (kin .eq. -1) then call nblank (basein, iflag2, n2) write (nouta,235) basein(1:n2)//'.xia', basein(1:n2)//'.eaa' 235 format (/, ' Files ''', a, ''' and ''', a, ''' are', + ' ASCII.') endif write (nouta,250) rident, pgm 250 format (/, ' From the header record(s):', + /, ' Data Identification =''', a56, '''', + /, ' Originating software =''', a8, '''') write (nouta,260) nc, nr 260 format (/, ' Number of columns =', i5, + ' Number of rows =', i5) write (nouta,270) ymin0, ymax0, dy 270 format (/, ' Latitude: minimum =', f8.0, + ' maximum =', f8.0, + ' increment =', f9.5) write (nouta,280) -xmax0, -xmin0, dx 280 format ( ' Longitude: minimum =', f8.0, + ' maximum =', f8.0, + ' increment =', f9.5) endif return end subroutine tonew (kin, kout, xmin, xmax, ymin, ymax, + ncfrst, nclast, nrfrst, nrlast, fouta, fouto) *** copy from input files to output files implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * nmax is the maximum size of the buffer. it is at least as large * as the largest number of columns in the input grid integer nmax parameter (nmax = 1024) * variables ending in an 'a' are associated with the xi grid * files while variables ending in an 'o' are associated with the * eta grid files. real xmin, xmax, ymin, ymax real angle, scr real bufa(nmax), bufo(nmax) integer ncfrst, nclast, nrfrst, nrlast integer nz, lrecl, kin, kout, nrpr1, nrpr2 integer iflag2, n2 integer i, i2, idum, irec1, irec2, j, jj, j1, j2, nr2, nc2 character*56 rident character*32 fouta, fouto character*8 pgm character*1 ascr integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto real xmin0, xmax0, ymin0, ymax0, dx, dy integer nc, nr common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc double precision zsuma, zsumo, zsqra, zsqro real zmina, zmino, zmaxa, zmaxo common /zstats/ zmina, zmaxa, zmino, zmaxo, + zsuma, zsqra, zsumo, zsqro data iflag2 /2/ * open output file, write header information, reopen for grid size angle = 0.0 rident = 'DEFLEC Extracted Region' pgm = 'DEFGRD' nz = 1 nr2 = nrlast - nrfrst + 1 nc2 = nclast - ncfrst + 1 if (kout .eq. 1) then * binary output files lrecl = 4*(nc2 + 1) open (nouta,file=fouta,form='unformatted',access='direct', + recl=lrecl,status='unknown') open (nouto,file=fouto,form='unformatted',access='direct', + recl=lrecl,status='unknown') write (nouta,rec=1) rident, pgm, nc2, nr2, nz, xmin, dx, + ymin, dy, angle write (nouto,rec=1) rident, pgm, nc2, nr2, nz, xmin, dx, + ymin, dy, angle elseif (kout .eq. -1) then * ascii transfer format output files open (nouta,file=fouta,form='formatted',access='sequential', + status='unknown') open (nouto,file=fouto,form='formatted',access='sequential', + status='unknown') write (nouta,140) rident, pgm write (nouto,140) rident, pgm 140 format (a56, a8) write (nouta,150) nc2, nr2, nz, xmin, dx, + ymin, dy, angle write (nouto,150) nc2, nr2, nz, xmin, dx, + ymin, dy, angle 150 format (3i4, 5f12.5) elseif (kout .eq. -2) then * ascii graphics format output files (surfer header information) open (nouta,file=fouta,form='formatted',access='sequential', + status='unknown') open (nouto,file=fouto,form='formatted',access='sequential', + status='unknown') write (nouta,160) write (nouto,160) 160 format ('DSAA') write (nouta,165) nc2, nr2 write (nouto,165) nc2, nr2 165 format (2i12) write (nouta,170) xmin, xmax write (nouto,170) xmin, xmax 170 format (2f12.4) write (nouta,170) ymin, ymax write (nouto,170) ymin, ymax write (nouta,170) zmina, zmaxa write (nouto,170) zmino, zmaxo endif call nblank (fouta, iflag2, n2) write (luout,80) fouta(1:n2), fouto(1:n2) 80 format (/, ' DEFGRD is now creating the ''', a, ''' and ''', a, + ''' files.', + /, ' This process may take several minutes.') * number of ascii rows per binary row (note that there are always * an odd number of grid columns) if (kin .eq. -1) then nrpr1 = nc/6 + 1 * for ascii input files, skip up to the nrfrst equivalent row rewind nina rewind nino read (nina, '(a1)') ascr read (nino, '(a1)') ascr read (nina, '(a1)') ascr read (nino, '(a1)') ascr if (nrfrst .gt. 1) then i2 = (nrfrst - 1)*nrpr1 do 210 i = 1, i2 read (nina, 200) scr read (nino, 200) scr 200 format (f12.6) 210 continue endif endif if (kout .eq. -1) then nrpr2 = (nclast - ncfrst + 1)/6 + 1 elseif (kout .eq. 1) then irec2 = 1 endif * copy input files to output files do 300 i = nrfrst, nrlast * read input file if (kin .eq. 1) then irec1 = i + 1 read (nina,rec=irec1) idum, (bufa(j), j = 1, nc) read (nino,rec=irec1) idum, (bufo(j), j = 1, nc) else do 320 jj = 1, nrpr1-1 j1 = (jj-1)*6 + 1 j2 = j1 + 5 read (nina, 310) (bufa(j), j = j1, j2) read (nino, 310) (bufo(j), j = j1, j2) 310 format (6f12.6) 320 continue j1 = j2 + 1 read (nina, 310) (bufa(j), j = j1, nc) read (nino, 310) (bufo(j), j = j1, nc) endif * write output file if (kout .eq. 1) then idum = 0 irec2 = irec2 + 1 write (nouta,rec=irec2) idum, (bufa(j), j = ncfrst, nclast) write (nouto,rec=irec2) idum, (bufo(j), j = ncfrst, nclast) elseif (kout .eq. -1) then do 420 jj = 1, nrpr2-1 j1 = (jj-1)*6 + ncfrst j2 = j1 + 5 write (nouta, 310) (bufa(j), j = j1, j2) write (nouto, 310) (bufo(j), j = j1, j2) 420 continue j1 = j2 + 1 write (nouta, 310) (bufa(j), j = j1, nclast) write (nouto, 310) (bufo(j), j = j1, nclast) elseif (kout .eq. -2) then write (nouta,440) (bufa(j), j = ncfrst, nclast) write (nouto,440) (bufo(j), j = ncfrst, nclast) 440 format (1024f12.6) endif 300 continue if (kout .eq. 1) then call nblank (fouta, iflag2, n2) write (luout,130) fouta(1:n2), fouto(1:n2) 130 format (/, ' Binary files ''', a, ''' and ''', a, '''', + ' have been created.') elseif (kout .eq. -1 .or. kout .eq. -2) then call nblank (fouta, iflag2, n2) write (luout,135) fouta(1:n2), fouto(1:n2) 135 format (/, ' ASCII files ''', a, ''' and ''', a, '''', + ' have been created.') endif return end subroutine zextrm (kin, kout, ncfrst, nclast, nrfrst, nrlast) * read from the input file to obtain the minimum and maximum z * values for the new grids, and other statistical values implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * nmax is the maximum size of the buffer. it is at least as large * as the largest number of columns in the input grid integer nmax parameter (nmax = 1024) * variables ending in an 'a' are associated with the xi grid * files while variables ending in an 'o' are associated with the * eta grid files. double precision avea, aveo, siga, sigo real bufa(nmax), bufo(nmax) real scr integer i2, nrpr1, jj, j1, j2, ntot integer irec, i, j, idum, kin, kout, ilast integer ncfrst, nclast, nrfrst, nrlast character*1 ascr integer luin, luout, nina, nino, nouta, nouto common /inout/ luin, luout, nina, nino, nouta, nouto double precision zsuma, zsumo, zsqra, zsqro real zmina, zmino, zmaxa, zmaxo common /zstats/ zmina, zmaxa, zmino, zmaxo, + zsuma, zsqra, zsumo, zsqro real xmin0, xmax0, ymin0, ymax0, dx, dy integer nc, nr common /grid0/ xmin0, xmax0, ymin0, ymax0, dx, dy, nr, nc * initialize zmina = 9.9e10 zmaxa = -9.9e10 zmino = 9.9e10 zmaxo = -9.9e10 zsuma = 0.d0 zsqra = 0.d0 zsumo = 0.d0 zsqro = 0.d0 ntot = 0 if (kout .ne. 0) then write (luout,80) 80 format (/, ' DEFGRD is now finding the maximum and minimum', + ' values in the new grids.',/, + ' This may take a few minutes.') else write (luout,85) 85 format (/, ' DEFGRD is now finding the maximum and minimum', + ' values in the new grids.',/, + ' This may take a few minutes.') endif if (kin .eq. 1) then do 300 i = nrfrst, nrlast irec = i + 1 read (nina,rec=irec) idum, (bufa(j), j = 1, nc) read (nino,rec=irec) idum, (bufo(j), j = 1, nc) call zsums (i, ncfrst, nclast, bufa, bufo) ntot = ntot + nclast - ncfrst + 1 300 continue elseif (kin .eq. -1) then nrpr1 = nc/6 + 1 * for ascii input files, skip up to the nrfrst equivalent row rewind nina rewind nino read (nina, '(a1)') ascr read (nino, '(a1)') ascr read (nina, '(a1)') ascr read (nino, '(a1)') ascr if (nrfrst .gt. 1) then i2 = (nrfrst - 1)*nrpr1 do 210 i = 1, i2 read (nina, 200) scr read (nino, 200) scr 200 format (f12.6) 210 continue endif ilast = nrlast - nrfrst + 1 do 340 i = 1, ilast do 320 jj = 1, nrpr1-1 j1 = (jj-1)*6 + 1 j2 = j1 + 5 read (nina, 310) (bufa(j), j = j1, j2) read (nino, 310) (bufo(j), j = j1, j2) 310 format (6f12.6) 320 continue j1 = j2 + 1 read (nina, 310) (bufa(j), j = j1, nc) read (nino, 310) (bufo(j), j = j1, nc) call zsums (i, ncfrst, nclast, bufa, bufo) ntot = ntot + nclast - ncfrst + 1 340 continue endif avea = zsuma / ntot aveo = zsumo / ntot siga = dsqrt( ( dble(ntot)*zsqra - zsuma*zsuma )/ + ( dble(ntot)*dble(ntot-1) ) ) sigo = dsqrt( ( dble(ntot)*zsqro - zsumo*zsumo )/ + ( dble(ntot)*dble(ntot-1) ) ) write (luout, 500) 500 format (/, 12x, 'DEFLECTION STATISTICS (seconds of arc)', + /, 1x, 71('-'), + /, 18x, (5x, 'Xi (N/S)', 5x, 'Eta (E/W)')) write (luout, 510) zmina, zmino 510 format (' Minimum', 14x, f9.3, 5x, f9.3) write (luout, 520) zmaxa, zmaxo 520 format (' Maximum', 14x, f9.3, 5x, f9.3) write (luout, 530) avea, aveo 530 format (/, ' Average', 14x, f9.3, 5x, f9.3) write (luout, 540) siga, sigo 540 format (' Standard Deviation', 3x, f9.3, 5x, f9.3) if (kout .eq. 0) then write (nouta, 500) write (nouta, 510) zmina, zmino write (nouta, 520) zmaxa, zmaxo write (nouta, 530) avea, aveo write (nouta, 540) siga, sigo write (luout,*) endif return end subroutine zsums (i, ncfrst, nclast, bufa, bufo) * do the summations for the subroutine zextrm implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) * nmax is the maximum size of the buffer. it is at least as large * as the largest number of columns in the input grid integer nmax parameter (nmax = 1024) * variables ending in an 'a' are associated with the xi grid * files while variables ending in an 'o' are associated with the * eta grid files. real bufa(nmax), bufo(nmax) real za, zo integer k, i, ncfrst, nclast double precision zsuma, zsumo, zsqra, zsqro real zmina, zmino, zmaxa, zmaxo common /zstats/ zmina, zmaxa, zmino, zmaxo, + zsuma, zsqra, zsumo, zsqro do 400 k = ncfrst, nclast *** seconds of arc statistics za = bufa(k) zo = bufo(k) zmina = min(zmina, za) zmaxa = max(zmaxa, za) zmino = min(zmino, zo) zmaxo = max(zmaxo, zo) zsuma = zsuma + dble(za) zsqra = zsqra + dble(za*za) zsumo = zsumo + dble(zo) zsqro = zsqro + dble(zo*zo) 400 continue return end real function rcard (chline, leng, ierr) *** read a real number from a line of card image. *** leng is the length of the card *** blanks are the delimiters of the real*4 variable implicit real (a-h, o-z) implicit integer (i-n) * implicit undefined (a-z) integer leng, ierr, i, j, ileng real var character*(*) chline ierr = 0 * find first non-blank character * do while line character is blank, i is first non-blank character i = 1 10 if ( chline(i:i) .eq. ' ' .or. chline(i:i) .eq. ',' ) then i = i + 1 * check for totally blank card if ( i .ge. leng) then rcard = 0.0e0 leng = 0 return endif goto 10 endif * find first blank character (or end of line) * do while line character is not a blank j = i + 1 20 if ( chline(j:j) .ne. ' ' .and. chline(j:j) .ne. ',' ) then j = j + 1 * check for totally filed card if ( j .gt. leng) then goto 40 endif goto 20 endif * j is now 1 more than the position of the last non-blank character 40 j = j - 1 * ileng is the length of the real string, it cannot be greater * than 15 characters ileng = j - i + 1 if (ileng .gt. 15) then stop 'rcard' endif * read the real variable from the line, and set the return var to it read (chline(i:j), 55, err=9999) var 55 format (f15.0) rcard = var * now reset the values of leng and chline to the rest of the card chline( 1 : leng ) = chline( (j+1) : leng ) leng = leng - j return * read error 9999 ierr = 1 return end