SUBROUTINE GRDCHK (POSX, POSY, INSIDE)

C
C ROUTINE CHECKS IF THE POINT HAVING COORDINATES (POSX, POSY)
C IS WITHIN THE REGION SPANNED BY THE GRID
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      LOGICAL INSIDE
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     1                ITREF
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCDGRD

      INSIDE = .TRUE.
      IF (NCDGRD .EQ. 0) RETURN

      IF (POSX .LT. GRDLX .OR. POSX .GT. GRDUX) THEN
         INSIDE = .FALSE.
      ENDIF
      IF (POSY .LT. GRDLY .OR. POSY .GT. GRDUY) THEN
         INSIDE = .FALSE.
      ENDIF

      RETURN
      END

CB::GRDPOS

      SUBROUTINE GRDPOS (POSX, POSY, I, J)

C
C********1*********2*********3*********4*********5*********6*********7**
C
C NAME:        GRDPOS
C VERSION:     9302.01   (YYMM.DD)
C WRITTEN BY:  MR. C. RANDOLPH PHILIPP
C PURPOSE:     THIS SUBROUTINE RETURNS THE COORDINATE OF THE LOWER
C              LEFT CORNER OF THE GRID CONTAINING POINT POSX, POSY 
C              
C  INPUT PARAMETERS FROM ARGUMENT LIST:
C  ------------------------------------
C POSX, POSY   THE POSITION IN THE GRIDED AREA 
C
C  OUTPUT PARAMETERS FROM ARGUMENT LIST:
C  -------------------------------------
C I, J         THE COORDINATES OF LOWER LEFT CORNER OF THE GRID
C              CONTAINING THE ABOVE POSITION
C
C  GLOBAL VARIABLES AND CONSTANTS:
C  -------------------------------
C NONE
C
C    THIS MODULE CALLED BY:   TCOORD
C
C    THIS MODULE CALLS:       NONE
C
C    INCLUDE FILES USED:      NONE
C
C    COMMON BLOCKS USED:      /CDGRID/
C
C    REFERENCES:  SEE RICHARD SNAY
C
C    COMMENTS:
C
C********1*********2*********3*********4*********5*********6*********7**
C    MOFICATION HISTORY:
C::9302.11, CRP, ORIGINAL CREATION FOR DYNAP
C********1*********2*********3*********4*********5*********6*********7**
CE::GRDPOS
     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF


      STEPX = (GRDUX - GRDLX) / ICNTX
      STEPY = (GRDUY - GRDLY) / ICNTY

      I = IDINT((POSX - GRDLX)/STEPX) + 1
      J = IDINT((POSY - GRDLY)/STEPY) + 1

      RETURN
      END



CB::GRDWEI

      SUBROUTINE GRDWEI (POSX, POSY, I, J, WEI)

C
C********1*********2*********3*********4*********5*********6*********7**
C
C NAME:        GRDWEI
C VERSION:     9302.01   (YYMM.DD)
C WRITTEN BY:  MR. C. RANDOLPH PHILIPP
C PURPOSE:     THIS SUBROUTINE RETURNS THE NORMALIZED WEIGHTS FOR 
C              A BI-LINEAR AVERAGE OVER A PLANE
C              
C  INPUT PARAMETERS FROM ARGUMENT LIST:
C  ------------------------------------
C POSX, POSY   THE POSITION IN THE GRIDED AREA 
C
C  OUTPUT PARAMETERS FROM ARGUMENT LIST:
C  -------------------------------------
C I, J         THE COORDINATES OF LOWER LEFT CORNER OF THE GRID
C              CONTAINING THE ABOVE POSITION
C WEI          A TWO BY TWO ARRAY CONTAINING THE NORMALIZED WEIGHTS
C              FOR THE CORNER VECTORS
C
C  GLOBAL VARIABLES AND CONSTANTS:
C  -------------------------------
C NONE
C
C    THIS MODULE CALLED BY:   TCOORD
C
C    THIS MODULE CALLS:       NONE
C
C    INCLUDE FILES USED:      NONE
C
C    COMMON BLOCKS USED:      /CDGRID/
C
C    REFERENCES:  SEE RICHARD SNAY
C
C    COMMENTS:
C
C********1*********2*********3*********4*********5*********6*********7**
C    MOFICATION HISTORY:
C::9302.11, CRP, ORIGINAL CREATION FOR DYNAP
C********1*********2*********3*********4*********5*********6*********7**
CE::GRDWEI
    
C**** COMPUTES THE WEIGHTS FOR AN ELEMENT IN A GRID

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION WEI(2,2)
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF	

C     COMPUTE THE LIMITS OF THE GRID CELL

      GRLX = GRDLX + (I - 1) * (GRDUX - GRDLX)/ICNTX
      GRUX = GRDLX + I * (GRDUX - GRDLX)/ICNTX
      GRLY = GRDLY + (J - 1) * (GRDUY - GRDLY)/ICNTY
      GRUY = GRDLY + J * (GRDUY - GRDLY)/ICNTY

C     CHECK THAT THE GIVEN POINT IS IN THE GRID CELL


      IF (.NOT.(GRLX.LE.POSX .AND. POSX.LE.GRUX)) THEN

C**** INSERT ERROR ROUTINE CALL HERE      

      ENDIF

      IF (.NOT.(GRLY.LE.POSY .AND. POSY.LE.GRUY)) THEN

C**** INSERT ERROR ROUTINE CALL HERE      

      ENDIF


C**** NORMALIZED WEIGHT FOR A POINT INSIDE THE GRID.

      DENOM = (GRUX - GRLX) * (GRUY - GRLY)

      WEI(1,1) = (GRUX - POSX) * (GRUY - POSY) / DENOM

      WEI(2,1) = (POSX - GRLX) * (GRUY - POSY) / DENOM

      WEI(1,2) = (GRUX - POSX) * (POSY - GRLY) / DENOM

      WEI(2,2) = (POSX - GRLX) * (POSY - GRLY) / DENOM

      RETURN
      END

CB::GRDVEC
C
      SUBROUTINE GRDVEC (I, J, VEL, B)
C
C********1*********2*********3*********4*********5*********6*********7**
C
C NAME:        GRDVEC
C VERSION:     9302.01   (YYMM.DD)
C WRITTEN BY:  MR. C. RANDOLPH PHILIPP
C PURPOSE:     THIS SUBROUTINE RETRIEVES THE APPROXIMATE VALUES OF THE
C              GRID NODE VELOCITIES FOR GRID (I,J) 
C              
C  INPUT PARAMETERS FROM ARGUMENT LIST:
C  ------------------------------------
C I, J         THE COORDINATES OF LOWER LEFT CORNER OF THE GRID
C              CONTAINING THE ABOVE POSITION
C B            THE ARRAY CONTAINING ALL THE APPROXIMATE VALUES
C              FOR THE ADJUSTMENT
C
C  OUTPUT PARAMETERS FROM ARGUMENT LIST:
C  -------------------------------------
C VEL          A TWO BY TWO ARRAY CONTAINING THE VELOCITY VECTORS
C              FOR THE CORNERS OF THE GRID
C
C  GLOBAL VARIABLES AND CONSTANTS:
C  -------------------------------
C NONE
C
C    THIS MODULE CALLED BY:   TCOORD
C
C    THIS MODULE CALLS:       NONE
C
C    INCLUDE FILES USED:      NONE
C
C    COMMON BLOCKS USED:      /CDGRID/
C
C    REFERENCES:  SEE RICHARD SNAY
C
C    COMMENTS:
C
C********1*********2*********3*********4*********5*********6*********7**
C    MOFICATION HISTORY:
C::9302.11, CRP, ORIGINAL CREATION FOR DYNAP
C********1*********2*********3*********4*********5*********6*********7**
CE::GRDVEC     


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION VEL(2,2,3), B(*)
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF	



      DO 30 II = 0,1
         DO 20 IJ = 0,1
            DO 10 IVEC = 1, 3
               INDEX = IUNGRD(I + II, J + IJ, IVEC)
               VEL(II + 1, IJ + 1, IVEC) = B(INDEX)
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE   

      RETURN
      END


CB::TCOORD
C
      SUBROUTINE TCOORD (GLAT0,GLON0,EHT0, B,
     &                   GLAT,GLON,EHT, ITIME, DTIME)
C
C********1*********2*********3*********4*********5*********6*********7**
C
C NAME:        TCOORD
C VERSION:     9302.01   (YYMM.DD)
C WRITTEN BY:  MR. C. RANDOLPH PHILIPP
C PURPOSE:     THIS SUBROUTINE COMPUTES A POSITION AT TIME ITIME GIVEN
C              THE APROXIMATE POSITION AT TIME ITREF AND THE VELOCITY
C              AT THE GIVEN POSITION
C              
C  INPUT PARAMETERS FROM ARGUMENT LIST:
C  ------------------------------------
C GLAT0,       THE POSITION AT EPOCH ITREF, DEFINED BY LATITUDE,
C GLON0,       LONGITUDE, AND HEIGHT
C GEHT0
C
C ITIME        THE TIME FOR THE NEW POSITION
C
C  OUTPUT PARAMETERS FROM ARGUMENT LIST:
C  -------------------------------------
C GLAT,       THE POSITION AT EPOCH ITIME, DEFINED BY LATITUDE,
C GLON,       LONGITUDE, AND HEIGHT
C GEHT
C
C  GLOBAL VARIABLES AND CONSTANTS:
C  -------------------------------
C NONE
C
C    THIS MODULE CALLED BY:   
C
C    THIS MODULE CALLS:       GRDPOS, GRDWEI, GRDVEC
C
C    INCLUDE FILES USED:      NONE
C
C    COMMON BLOCKS USED:      /CDGRID/
C
C    REFERENCES:  SEE RICHARD SNAY
C
C    COMMENTS:
C
C********1*********2*********3*********4*********5*********6*********7**
C    MOFICATION HISTORY:
C::9302.11, CRP, ORIGINAL CREATION FOR DYNAP
C********1*********2*********3*********4*********5*********6*********7**
CE::TCOORD
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION WEI(2,2), VEL(2,2,3)
      DIMENSION B(*)
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF	
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

      IF (NCD.GT.0) THEN
C**** GET THE COORDINATES OF THE LOWER LEFT HAND CORNER OF THE GRID

        CALL GRDPOS(GLON0, GLAT0, I, J)

C**** GET THE WEIGHTS FOR THE FOUR CORNERS

        CALL GRDWEI(GLON0, GLAT0, I, J, WEI)

C**** GET THE VELOCITY VECTORS AT THE FOUR CORNERS

        CALL GRDVEC(I, J, VEL, B)

C**** COMPUTE THE TIME DIFFERENCE BETWEEN THE REFERENCE EPOCH
C**** AND THE ITIME

        DTIME = DBLE(ITIME - ITREF) / (365.25D0 * 24.0D0 * 60.0D0)

C**** COMPUTE THE NEW POSITION AT TIME ITIME

C****
C**** THE FOLLOWING EQUATIONS COMPUTES THE VELOCITY AT POSITION
C**** (GLAT0, GLON0, EHT0).  USING THE NORMALIZED WEIGHTS FROM
C**** FUNCTION GRDWEI, AN AVERAGE VELOCITY IS COMPUTED.  USING
C**** THE AVERAGE VELOCITY, THE NEW POSITION IS COMPUTED FOR TIME
C**** ITIME.
C****

        GLAT = GLAT0 + (WEI(1,1) * VEL(1,1,1)
     &               +  WEI(1,2) * VEL(1,2,1)
     &               +  WEI(2,1) * VEL(2,1,1)
     &               +  WEI(2,2) * VEL(2,2,1)) * DTIME
      
        GLON = GLON0 + (WEI(1,1) * VEL(1,1,2)
     &               +  WEI(1,2) * VEL(1,2,2)
     &               +  WEI(2,1) * VEL(2,1,2)
     &               +  WEI(2,2) * VEL(2,2,2)) * DTIME
 
        EHT = EHT0   + (WEI(1,1) * VEL(1,1,3)
     &               +  WEI(1,2) * VEL(1,2,3)
     &               +  WEI(2,1) * VEL(2,1,3)
     &               +  WEI(2,2) * VEL(2,2,3)) * DTIME
      ELSE
        EHT = EHT0
        GLON = GLON0
        GLAT = GLAT0
        DTIME = 0.D0
      ENDIF
 
      RETURN
      END

CB::VELOC
C
      SUBROUTINE VELOC (GLAT,GLON,B,VN,VE,VU,SN,SE,CORR,A,NX,SIGUWT)
C
C****************************************************************
C
C NAME:         VELOC
C VERSION:      9310.06
C WRITTEN BY:   RICHARD SNAY
C PURPOSE:      THIS SUBROTINE COMPUTES THE ERROR ELLIPSE
C               AT A GIVEN LOCATION
C
C  INPUT PARAMETERS FROM ARGUMENT LIST:
C  ____________________________________
C  GLAT         LATITUDE OF POSITION AT REFERENCE EPOCH
C  GLON         LONGITUDE OF POSITION AT REFERENCE EPOCH
C  B            ARRAY CONTAINING ESTIMATES OF GRID VELOCITIES
C
C  OUTPUT PARAMETERS FROM ARGUMENT LIST:
C  --------------------------------------
C  VN           NORTHWARD VELOCITY in mm/yr
C  VE           EASTWARD VELOCITY in mm/yr
C  VU           UPWARD VELOCITY in mm/yr
C  SN		STANDARD ERROR FOR NORTHWARD
C		VELOCITY in mm/yr
C  SE		STANDARD ERROR FOR EASTWARD
C 		VELOCITY in mm/yr
C  CORR		CORRELATION BETWEEN
C		NORTHWARD AND EASTWARD
C		VELOCITIES (unitless)
C
C  GLOBAL VARIABLES AND CONSTANTS:
C  -------------------------------------
C  NONE
C
C    THIS MODULE IS CALLED BY:    ADJPOS
C
C    THIS MODULE CALLS:       GRDPOS, GRDWEI, GRDVEC, RADCUR
C
C    INCLUDE FILES USED:          NONE
C
C    COMMON BLOCK USED:           /CDGRID/
C
C    REFERENCES:    SEE RICHARD SNAY
C
C    COMMENTS:
C
C************************************************************
C    MODIFICATION HISTORY: C::931006, RAS, ORIGINAL CREATION FOR DYNAPG 
C************************************************************
CE::VELOC
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION WEI(2,2), VEL(2,2,3)
      DIMENSION B(*),A(*),NX(*)
      DIMENSION II(4),JJ(4),COEF(1,4)
      LOGICAL LMSL,LSS,LUP
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /CONST/ PI,PI2,RAD,RADSEC,TWOPI

C**** GET THE COORDINATES OF THE LOWER LEFT CORNER OF THE GRID

      CALL GRDPOS (GLON,GLAT, I, J)

C**** GET THE WEIGHTS FOR THE FOUR CORNERS

      CALL GRDWEI (GLON, GLAT, I, J, WEI)

C**** GET THE VELOCITY VECTORS AT THE FOUR CORNERS

      CALL GRDVEC (I, J, VEL, B)

C**** COMPUTE VELOCITY VECTOR AT LOCATION

      VN = WEI(1,1) * VEL(1,1,1) + WEI(1,2) * VEL(1,2,1)
     &   + WEI(2,1) * VEL(2,1,1) + WEI(2,2) * VEL(2,2,1)

      VE = WEI(1,1) * VEL(1,1,2) + WEI(1,2) * VEL(1,2,2)
     &   + WEI(2,1) * VEL(2,1,2) + WEI(2,2) * VEL(2,2,2)

      VU = WEI(1,1) * VEL(1,1,3) + WEI(1,2) * VEL(1,2,3)
     &   + WEI(2,1) * VEL(2,1,3) + WEI(2,2) * VEL(2,2,3)

C**** CONVERT VELOCITIES TO MM/YR

      CALL RADCUR (GLAT,RMER,RPV)
      RPAR = RPV * DCOS(GLAT)
      VN = VN * RMER * 1000.D0
      VE = VE * RPAR * 1000.D0
      VU = VU * 1000.D0

C**** COMPUTING ERROR ELLIPSE

      IF (IMODE .NE. 3) THEN
        SN = 0.0D0
        SE = 0.0D0
        CORR = 0.0D0
      ELSE
	II(1)=I
	II(2)=I+1
	II(3)=I
	II(4)=I+1
	JJ(1)=J
	JJ(2)=J
	JJ(3)=J+1
	JJ(4)=J+1
	COEF(1,1)=WEI(1,1)
	COEF(1,2)=WEI(2,1)
	COEF(1,3)=WEI(1,2)
	COEF(1,4)=WEI(2,2)

	KVEC=1
	LVEC=1
	CALL GETCOV(II,JJ,KVEC,LVEC,A,NX,COEF,COV)
	SN=DSQRT(COV)

	KVEC=2
	LVEC=2
	CALL GETCOV(II,JJ,KVEC,LVEC,A,NX,COEF,COV)
	SE=DSQRT(COV)

	KVEC=1
	LVEC=2
	CALL GETCOV(II,JJ,KVEC,LVEC,A,NX,COEF,COV)

	CORR=COV/(SN*SE)
	SN=SN*RMER*1000.D0/RADSEC
	SE=SE*RPAR*1000.D0/RADSEC
	IF (LSS) THEN
	  SN=SN*SIGUWT
	  SE=SE*SIGUWT
	ENDIF
      ENDIF 
      RETURN
      END

      SUBROUTINE GETCOV(II,JJ,KVEC,LVEC,A,NX,COEF,COV)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      LOGICAL GETA
      DIMENSION A(*),NX(*)
      DIMENSION II(4),JJ(4),COEF(1,4)
      DIMENSION ACOV(4,4),RCOV(1,1),W(4)

      DO K=1,4
	INDEXK=IUNGRD(II(K),JJ(K),KVEC)
	DO L=1,4
	  INDEXL=IUNGRD(II(L),JJ(L),LVEC)
	  IF (.NOT.GETA(INDEXK,INDEXL,ACOV(K,L),A,NX)) THEN
	    WRITE (6,100)
  100	    FORMAT ('0ERRROR IN GETTING VELOCITY COVARIANCE') 
	    STOP
	  ENDIF
	ENDDO
      ENDDO

      CALL ABAT(COEF,ACOV,RCOV,W,1,4)
      COV=RCOV(1,1)
      RETURN
      END

      SUBROUTINE AFGR (CARD,B,FATAL)
*** READ CRUSTAL MOTION REFERENCE FRAME PARAMETERS

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      CHARACTER*80 CARD, CARDGR
      CHARACTER*1 TCREF
      CHARACTER*10 LONMIN,LONMAX,LATMIN,LATMAX
      LOGICAL FATAL
      DIMENSION B(*)    

      COMMON /CARDG/ CARDGR
      COMMON /CONST/ PI,PI2,RAD,RADSEC,TWOPI
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF

      IF (NAUX.GT.0) THEN
          WRITE(6,666)
  666     FORMAT('0ERROR - GR CARD MUST PRECEED ALL SS CARDS.',/)
          CALL ABORT2
      ENDIF

      CARDGR = CARD

      READ (CARD,10) IYRREF,IMOREF,IDYREF,
     &              LONMIN, LONMAX, ICNTX,
     &              LATMIN, LATMAX, ICNTY

   10 FORMAT(2X,I4, I2, I2, 2(2A10, I3))  

      CALL RDEG(LONMIN, GRDLX, 'W')
      CALL RDEG(LONMAX, GRDUX, 'W')
      CALL RDEG(LATMIN, GRDLY, 'S')
      CALL RDEG(LATMAX, GRDUY, 'S')

      IHRREF = 0
      IMNREF = 0
      TCREF = 'Z'
      
*** SET DEFAULTS

      CALL TOMNT (IYRREF,IMOREF,IDYREF,IHRREF,IMNREF,TCREF,ITREF)

      NCD = 3 * (ICNTX + 1) * (ICNTY + 1)

      DO 20 I = 1, NCD
         B(I) = 0.0D0
   20 CONTINUE

      RETURN
      END

      SUBROUTINE RDEG (INPUT,VAL,CNEG)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      CHARACTER*10 INPUT
      CHARACTER*1  CNEG
      INTEGER      DEG, MIN, ISEC
      COMMON /CONST/ PI, PI2, RAD, RADSEC, TWOPI

      DO 10 I = 1, 9
         IF (INPUT(I:I).EQ.' ') THEN
             INPUT(I:I) = '0'
         ENDIF
   10 CONTINUE

      READ (INPUT,20) DEG, MIN, ISEC

   20 FORMAT(I3,I2,I4)

      SEC = ISEC/100.D0

      VAL = (DEG + (MIN/60.D0) + (SEC/3600.D0) ) / RAD

      IF (INPUT(10:10).EQ.CNEG) THEN
         VAL = -VAL
      ENDIF

      RETURN
      END



      INTEGER FUNCTION IUNGRD(I, J, IVEC)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCDGRD
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF

      IUNGRD = 3 * ((J - 1) * (ICNTX + 1) +  (I - 1)) + IVEC

      RETURN
      END
     
      SUBROUTINE GRDLOC (I, J, POSX, POSY)

*** DETERMINE LOCATION OF I,J-GRID NODE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY,
     &                ITREF

      STEPX = (GRDUX - GRDLX) / ICNTX
      STEPY = (GRDUY - GRDLY) / ICNTY

      POSX = GRDLX + (I - 1)*STEPX
      POSY = GRDLY + (J - 1)*STEPY

      RETURN
      END