      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 
