invwt.f

Go to the documentation of this file.

C*********************************************************************
    SUBROUTINE rinvmod(TRPVAL,X0,Y0,XX,YY,NOBS,VARBL)
C*********************************************************************
C 
C
C  MODIFIED INVERSE DISTANCE INTERPOLATION;  USES ONLY NEARBY
C  OBSERVATIONS.  THE WEIGHTING FUNCTION IS WNDWT (X0,Y0,X1,Y1).
C  INTERPOLATES TO THE POINT X0,Y0 FROM VALUES OF THE VARIABLE (VARBL)
C  AT POINTS XX, YY.  IF A SINGLE OBSERVED VALUE IS WITHIN THE MINIMUM
C  DISTANCE SPECIFIED IN WNDWT, IT IS USED.  ASSUMES THAT GRID UNITS 
C  ARE USED.
C
C       FLUDWIG, 5/2002
C
        include 'NGRIDS.PAR'
C
        include 'LIMITS.CMM'
C
    parameter(npts=26)
C
    REAL VARBL(nsites),XX(nsites),YY(nsites)
C
    d2min=0.1
C
C  IF ONLY ONE OBS, USE THAT VALUE EVERYWHERE.
C
    IF (nobs.EQ.1) THEN
       trpval=varbl(1)
       RETURN
    ELSE IF (nobs.EQ.0) THEN
      WRITE (*,*) 'SOMETHING WRONG -- NOBS= ',nobs
      stop
    END IF
C
C  CHECK FOR TOO MANY INPUT VALUES
C
    IF (nobs.GT.nsites) THEN
       WRITE (*,*) 'TOO MANY OBS IN RM2MOD; USING 1ST ',nsites
       ntry=nsites
    ELSE
       ntry=nobs
    END IF
C
C  DEFINE EVER INCREASING RADIUS UNTIL ENOUGH OBS AVAILABLE FOR
C  INTERPOLATION.  MAXIMUM SET TO BE TWICW THE MAXIMUM RADIUS OF THE 
C  DOMAIN (RMAX).  SEARCHED AREA INCREASES BY RSTEP UNTIL AT LEAST 2
C  OBSERVATIONS FOUND.  RSTEP IS 1/20 MAX NUMBER OF GRID UNITS.
C  FIRST STEP IS A RADIUS OF SQRT(2) GRID UNITS. IF ANY VALUE FOUND 
C  THIS  CLOSE, IT IS USED.
C
    rmax=sqrt(float(nxgrd**2+nygrd**2))
    rstep=rmax/20.0
C
C  3 CHOICES EXAMINED FOR RMIN:
C       1.  RADIUS (0.5) OF CIRCLE TANGENT TO THE 4 SIDES OF A 
C           GRID SQUARE
C       2.  RADIUS (0.707) OF CIRCLE THROUGH THE 4 CORNERS OF A 
C           GRID SQUARE
C   3.  RADIUS (.604) MIDWAY BETWEEN 0.5 & SQRT(1/2)
C       4.  IN SQUARE OF 0.5 UNIT SIDE WITH CORNER AT GRID SQUARE
C
C  APPROACH 3 SEEMS TO GET MOST OF OBSERVATIONS, WHILE AVERAGING 
C  MORE THAN ONE OB ON FEWER OCCASIONS.
C
    rmin=0.604
C
    sum=0.0
    sumwt=0.0
    ncase=0
    r2past=0.0
    DO 200 nr=0,20
       rad=rmin+float(nr)*rstep
       r2test=rad**2
       DO 150 iob=1,ntry
          IF (nint(varbl(iob)).NE.-9999) THEN
             r2=(x0-xx(iob))**2+(y0-yy(iob))**2
C
C  IF OBSERVATION IS WITHIN A GRID SQUARE AROUND THIS POINT, USE IT.
C  OTHERWISE, GET VALUES FOR INVERSE DISTANCE INTRPOLATION.
C
             IF (nr.EQ.0 .AND. r2.LE.r2test) THEN
C
            wt=wndwt(x0,y0,xx(iob),yy(iob))
            sum=sum+wt*varbl(iob)
            sumwt=sumwt+wt
             ELSE IF (r2.LE.r2test .AND. r2.GT.r2past) THEN
C
            wt=wndwt(x0,y0,xx(iob),yy(iob))
            sum=sum+wt*varbl(iob)
            sumwt=sumwt+wt
            ncase=ncase+1
         END IF
          END IF
150    CONTINUE
C
C  TE FOLLOWING USES NEARBY OBS WITH LITTLE OR NO INTERPOLATION.
C
       IF (nr.EQ.0 .AND. sumwt .GT. 0.0) THEN
          trpval=sum/sumwt
          RETURN
C
C  WHEN WE HAVE ACCUMULATED AT LEAST 15 PTS, STOP ADDING POINTS 
C  AND INTERPOLATE 
C
       ELSE IF (ncase .GE. 15) THEN
C
          trpval=sum/sumwt
          RETURN
       END IF
       r2past=r2test
200 CONTINUE
C
C  WE NEVER GOT 15 PTS, SO WE USE THE INVERSE DISTANCE WEIGHTED 
C  ESTIMATE WITH WHAT WE HAVE.
C
    IF (sumwt.GT.0.0) THEN
       trpval=sum/sumwt
       RETURN
    ELSE
       WRITE (*,*) 'SOMETHING WRONG -- NCASE= ', ncase
       stop
    END IF
C
    END
C
C*********************************************************************
    SUBROUTINE rinvmold(TRPVAL,X0,Y0,XX,YY,NOBS,VARBL)
C*********************************************************************
C                  
C  INVERSE DISTANCE INTERPOLATION;  USES ALL AVAILABLE
C  OBSERVATIONS.  THE WEIGHTING FUNCTION IS WNDWT (X0,Y0,X1,Y1).
C  INTERPOLATES TO THE POINT X0,Y0 FROM VALUES OF THE VARIABLE (VARBL)
C  AT POINTS XX, YY.  IF A SINGLE OBSERVED VALUE IS WITHIN THE MINIMUM
C  DISTANCE SPECIFIED IN WNDWT, THE FIELD AROUND THAT PT IS ADJUSTED 
C  ELSEWHERE TO MATCH THE VALUE AT THE PT.
C
C       FLUDWIG, 7/2002
C
        include 'NGRIDS.PAR'
C
        include 'LIMITS.CMM'
C
    parameter(npts=nsites)
C
    REAL VARBL(nsites),XX(nsites),YY(nsites)
C
C  IF ONLY ONE OBS, USE THAT VALUE EVERYWHERE.
C
    IF (nobs.GT.1 .AND. nobs.LE. nsites) THEN
       sum=0.0
       sumwt=0.0
       ncase=0
       DO 200 iob=1,nobs
          wt=wndwt(x0,y0,xx(iob),yy(iob))
          sum=sum+wt*varbl(iob)
          sumwt=sumwt+wt
200    CONTINUE
       trpval=sum/sumwt
       RETURN
    ELSE IF (nobs.EQ.1) THEN
       trpval=varbl(1)
       RETURN
    ELSE IF (nobs.GT.nsites) THEN
       WRITE (*,*) 'TOO MANY OBS IN RM2MOLD '
       stop
    ELSE IF (nobs.EQ.0) THEN
      WRITE (*,*) 'SOMETHING WRONG IN RM2MOLD -- NOBS= ',nobs
      stop
    END IF
C
    END