program vertcon *********************************************************************** * Release 2.0 * * PURPOSE: CALCULATE THE DATUM CONVERSION CORRECTION TO NGVD29 * * ORTHOMETRIC HEIGHT IN ORDER TO OBTAIN NAVD88 HEIGHT. * * THE INPUT REQUIRES LATITUDE AND LONGITUDE VALUES; * * * * [THE GEODETIC DATUM TO WHICH THE POSITION IS REFERENCED * * (I.E., NAD27 OR NAD83) IS GENERALLY IRRELEVANT DUE TO * * THE IMPRECISION OF SCALED NAD27 POSITIONS; HOWEVER, * * SINCE NUMEROUS BENCH MARKS ARE ALSO HORIZONTAL NETWORK * * STATIONS, ALL POSITIONS WILL BE CONVERTED TO NAD83 * * POSITIONS IN THE NEXT RELEASE] * * * * THE COMPUTATION IS PERFORMED AS AN INTERPOLATION USING * * A BILINEAR INTERPOLATOR AMONG FOUR GRID POINTS. * * * * THE PROGRAM REQUIRES THAT THE USER SPECIFY: * * 1) THE NAME OF AN INPUT FILE (IF AVAILABLE). * * 2) THE NAME OF AN OUTPUT FILE * * * * THIS PROGRAM ALLOWS FOR : * * GENERIC FILE FORMATS - SEE SUBROUTINE TYPE1 OR TYPE2 * * OR * * NGS INTERNAL BENCH MARK FILE FORMAT - SEE TYPE3 * * * * [THE CODE CAN BE EASILY MODIFIED TO ACCOMMODATE CUSTOM * * FILE SPECIFICATIONS (TYPE..) BY MODIFYING SUBROUTINES: * * GETPT, IPARMS, WRTPT, AND (OPTIONALLY) FHELP.] * * * * * * VERSION CODE: 2.00 (August-94) * * * * INFORMATION CONCERNING THE 'VERTCON' SYSTEM MAY BE OBTAINED FROM: * * RUDOLF J. FURY NOAA/NOS/C&GS/NGS * * SILVER SPRING, MARYLAND 20910-3282 * *********************************************************************** * * * 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. * * * *********************************************************************** * c real vrsion parameter (vrsion = 2.00E0) logical page, screen,nogo character*60 name SAVE ITYPE,PAGE,SCREEN * VARIABLES FOR THE CORPSCON CONVERSION * CORPSD - ARRAY HOLDING ALL POINT INFORMATION * FLAG - FLAG FOR SUCCESS/FAILURE OF CONVERSION * INFILE - INPUT FILE NAME WITH CORPSD DATA * NAME - POINT NAME DOUBLE PRECISION CORPSD(29) INTEGER*2 KEY,FLAG,HDATUM,IUNIT,OUNIT,METER CHARACTER*60 INFILE CHARACTER*25 PTNAME * VARIABLES FROM MLOOP ROUTINE. MLOOP IS SKIPPED. DOUBLE PRECISION XPT,YPT SAVE IOS REAL GHT REAL*8 DHPRED,XEAST c COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, LDUMP, NSPACE(2) * INITIALIZE VARIABLES CALL INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE) * LOAD 'VERTCON' GRID WRITE(*,*) 'CORPSVRT STARTED' call loadgrd(nogo) if(.not.nogo) then READ(*,10) INFILE 10 FORMAT(A60) READ(*,15) HDATUM READ(*,15) KEY READ(*,15) IUNIT READ(*,15) OUNIT 15 FORMAT(I1) METER = 3 WRITE(*,*) 'CORPSVRT INPUT PARAMETERS READ SUCCESSFULLY' OPEN(81,FILE=INFILE,ACCESS='SEQUENTIAL',FORM='BINARY', & ERR=25,STATUS='OLD') WRITE(*,*) 'CORPSVRT INPUT FILE OPEN SUCCESSFULLY' OPEN(83,FILE='INT$OUT',ACCESS='SEQUENTIAL',FORM='BINARY', & ERR=25,STATUS='UNKNOWN') WRITE(*,*) 'CORPSVRT OUTPUT FILE OPEN SUCCESSFULLY' 20 READ(81,END=25) PTNAME,(CORPSD(II),II=1,29),FLAG C HDATUM=0=NAD83, HDATUM=1=NAD27 IF(HDATUM.EQ.0) THEN YPT = CORPSD(11) XPT = CORPSD(12) ELSE YPT = CORPSD(9) XPT = CORPSD(10) ENDIF XEAST = 360.0-XPT CALL INTERP(YPT,XEAST,DHPRED,IOS) IF(IOS.NE.0) FLAG = FLAG + 16 GHT=DHPRED*0.001D0 GHT=ROUND(GHT) C KEY=0=29TO88, KEY=1=88TO29 C FOR 29TO88 ADD GHT, FOR 88TO29 SUB GHT IF(KEY.EQ.0) THEN IF(IUNIT.NE.3) CALL CVT(IUNIT, METER, CORPSD(23)) CORPSD(24) = CORPSD(23) + GHT IF(IUNIT.NE.3) CALL CVT(METER, IUNIT, CORPSD(23)) IF(OUNIT.NE.3) CALL CVT(METER, OUNIT, CORPSD(24)) ELSE IF(IUNIT.NE.3) CALL CVT(IUNIT, METER, CORPSD(24)) CORPSD(23) = CORPSD(24) - GHT IF(IUNIT.NE.3) CALL CVT(METER, IUNIT, CORPSD(24)) IF(OUNIT.NE.3) CALL CVT(METER, OUNIT, CORPSD(23)) ENDIF WRITE(83) PTNAME,(CORPSD(II),II=1,29),FLAG GOTO 20 endif 25 CLOSE(81) CLOSE(83) WRITE(*,*) 'CORPSVRT INPUT/OUTPUT FILES CLOSED' WRITE(*,*) 'CORPSVRT COMPLETE' STOP C write(LUOUT,*) 'Normal End' stop 1001 write(LUOUT,1104) ios 1104 format(' Input File i/o err. =',i5) stop 1002 write(LUOUT,1105) ios 1105 format(' Output File i/o err. =',i5) stop 1003 write(LUOUT,1106) ios 1106 format(' Dump File i/o err. =',i5) stop end * * THIS SUBROUTINE CONVERTS BETWEEN U.S. SURVEY FEET, INTERNATIONAL * FEET, AND METERS. CURR IS INPUT, DES IS OUTPUT, AN PT IS POINT * TO CONVERT. US FOOT-1, INTL FOOT-2, METER-3. THIS IS TAKEN FROM * DIRECTLY FROM CORPSCON SOURCE CODE. * SUBROUTINE CVT(CURR, DES, PT) DOUBLE PRECISION PT DOUBLE PRECISION USF, USM, INF, INM INTEGER*2 CURR, DES IF(CURR.EQ.DES) RETURN USF = 3937.0 USM = 1200.0 INF = 1250.0 INM = 381.0 IF(CURR.EQ.1) THEN PT = PT * USM / USF IF(DES.EQ.2) PT = PT * INF / INM RETURN ENDIF IF(CURR.EQ.2) THEN PT = PT * INM / INF IF(DES.EQ.1) PT = PT * USF /USM RETURN ENDIF IF(CURR.EQ.3) THEN IF(DES.EQ.1) THEN PT = PT * USF / USM ELSE PT = PT * INF / INM ENDIF ENDIF RETURN END SUBROUTINE INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE) parameter (MAXAREA=10,NSPACE1=6*MAXAREA,NSPACE2=2*MAXAREA) CHARACTER*12 FILES CHARACTER*20 B20 CHARACTER*80 B80 CHARACTER*80 CARD CHARACTER*40 NAME CHARACTER*96 B96 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) c INTEGER IPAGE, ITYPE LOGICAL PAGE, SCREEN REAL*4 SPACE1,SPACE2,MARGIN COMMON /CURNT/ B96 COMMON /GRIDS/ FILES(MAXAREA) COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, LDUMP, NSPACE(2) common/GSTUFF/ SPACE1(NSPACE1),MARGIN(MAXAREA),SPACE2(NSPACE2) EQUIVALENCE(CARD,B96) * Initialize card variable in common CURNT to blank CARD = B80 * Set the logical units for input/output common INOUT LUIN = 5 LUOUT = 6 NOUT = 101 NIN = 102 LDUMP = LUOUT * 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 FILES(1)='vertconw.94 ' MARGIN(1)=5.0 FILES(2)='vertconc.94 ' MARGIN(2)=5.0 FILES(3)='vertcone.94 ' MARGIN(3)=0.0 DO N=4,MAXAREA FILES(N)=' ' ENDDO RETURN END REAL FUNCTION ROUND(VALUE) REAL VALUE INTEGER INTHT INTHT=(SIGN(1.0,VALUE)*(ABS(VALUE)+0.0005))*1000. ROUND=INTHT*0.001 RETURN END subroutine angle(degrees,deg,min,sec) real*8 degrees integer*4 deg deg=degrees min=(degrees-deg)*60.d0 sec=(degrees-deg-min/60.d0)*3600.d0 if((sec+0.000001).ge.60.0) then min=min+1 sec=sec-60.0 endif if(min.ge.60) then deg=deg+1 min=min-1 endif return end subroutine interp(gla,glo,val,ios) *** bilinear interpolation from a standard grid file implicit double precision(a-h,o-z) parameter(narea=10,len=461) logical nogo integer nz(1) real*4 glamn,dgla,glamx,glomn,dglo,glomx,margin,skip character*8 pgm character*12 files character*56 ident common/grids/ files(narea) common/gstuff/glamn(narea),dgla(narea),glamx(narea),glomn(narea), . dglo(narea),glomx(narea),margin(narea),nla(narea),nlo(narea) *** do the interpolation, check for clipped points ios=0 lin =90 do n=1,narea lin=lin+1 c if(gla.ge.glamn(n).and.gla.le.glamx(n).and. * glo.ge.glomn(n).and.glo.le.(glomx(n)-margin(n))) then c call fgrid(lin,glo,gla,glomn(n),glamn(n), . dglo(n),dgla(n),nlo(n),xgrid,ygrid,irow,jcol, . tee1,tee2,tee3,tee4,ios) if(ios.ne.0) return c call coeff(tee1,tee2,tee3,tee4,ay,bee,cee,dee) call surf(xgrid,ygrid,val,ay,bee,cee,dee,irow,jcol) c return endif enddo c ios=999 return entry loadgrd(nogo) lin =90 nogo=.true. c open available grid files do n=1,narea if(files(n).ne.' ') then lin=lin+1 nogo=.false. open(lin,file=files(n),status='old',form='unformatted', . access='direct',recl=1848) WRITE(*,*) files(n),' OPEN SUCCESSFULLY' read(lin,rec=1) ident,pgm,nlo(n),nla(n),nz(n), . glomn(n),dglo(n),glamn(n),dgla(n),skip c write(6,'(i4,4x,a12/1x,a56,a8/3i5,5f10.5)') c . n,files(n),ident,pgm,nlo(n),nla(n),nz(n), c . glomn(n),dglo(n),glamn(n),dgla(n),margin(n) c WRITE (6,933) c 933 FORMAT(20x,'(Hit RETURN to continue)') c READ (5,'(A1)') ANS c glamx(n)=glamn(n)+dgla(n)*(nla(n)-1) glomn(n)=360.0+glomn(n) glomx(n)=glomn(n)+dglo(n)*(nlo(n)-1) if(nlo(n).gt.len) stop 12345 if(nla(n).lt.3.or.nlo(n).lt.3) stop 54321 else margin(n)=0.0 glamn(n)=0.0 glomn(n)=0.0 glamx(n)=0.0 glomx(n)=0.0 endif enddo return entry closegrd lin=90 do n=1,narea if(files(n).ne.' ') then lin=lin+1 close(lin,status='keep') WRITE(*,*) 'CORPSVRT DATA FILES CLOSED' endif enddo return end subroutine fgrid(lin,glo,gla,glomn,glamn,dglo,dgla,nlo, . xgrid,ygrid,irow,jcol,tee1,tee2,tee3,tee4,ios) implicit double precision (a-h,o-z) parameter(len=461) integer dummy real*4 z(len) real*4 glamn,dgla,glomn,dglo c calculate the coordinates for the point in terms of grid c indices xgrid=(glo-glomn)/dglo+1.d0 ygrid=(gla-glamn)/dgla+1.d0 c find the i,j, values for the SW corner of local square irow=int(ygrid) jcol=int(xgrid) c read(lin,rec=irow+1) dummy,(z(i),i=1,nlo) if((z(jcol).eq.9999.0).or.(z(jcol+1).eq.9999.0)) go to 10 c tee1=z(jcol) tee3=z(jcol+1) c read(lin,rec=irow+2) dummy,(z(i),i=1,nlo) if((z(jcol).eq.9999.0).or.(z(jcol+1).eq.9999.0)) go to 10 tee2=z(jcol) tee4=z(jcol+1) return c 10 ios=999 return end subroutine coeff (tee1,tee2,tee3,tee4,ay,bee,cee,dee) implicit double precision (a-h,o-z) ay=tee1 bee=tee3-tee1 cee=tee2-tee1 dee=tee4-tee3-tee2+tee1 return end c c********************************************************************* c* subroutine surf: interpolates the z value * c********************************************************************* c subroutine surf (xgrid,ygrid,zee,ay,bee,cee,dee,irow,jcol) c c calculates the value of the grid at the point xpt, ypt. the interp. c is done in the index coordinate system for convenience. implicit double precision (a-h,o-z) implicit integer (i-n) c zee1 = ay zee2 = bee*(xgrid - float(jcol)) zee3 = cee*(ygrid - float(irow)) zee4 = dee*(xgrid - float(jcol))*(ygrid - float(irow)) zee = zee1+zee2+zee3+zee4 return end