interp.f

Go to the documentation of this file.

C********************************************************************
        SUBROUTINE sfctrp(NCALL)
C********************************************************************
C
C  THIS SUBROUTINE GETS FIRST ESTIMATE OF SELECTED MET VARIABLES
C  AT LOWEST GRID POINTS. USES INVERSE DISTANCE-SQUARED WHEN FALSE.
C
C   THIS VERSION FLUDWIG 3/2002
C
        include 'NGRIDS.PAR'
C
        include 'ANCHOR.CMM'
        include 'FLOWER.CMM'
        include 'LIMITS.CMM'
        include 'STALOC.CMM'
C
    REAL UOB(nsites),VOB(nsites),XVAL(nsites),YVAL(nsites)
C
    SAVE
C
        DATA zero /0.0/
C
C  FIRST, GET A VALUE FOR EACH X,Y GRID POINT.
C
C  INTERPOLATE TO GET VALUES AT 10 M U(IX,IY,1) & V(IX,IY,1)
C  ABOVE EACH GRID POINT ON 1ST CALL
C
    IF (ncall .LE.1) THEN
C
C*********BUG FIX WHEN MULTIPLE CASES RUN*****1/2002********
C  FLOW SURFACE HEIGHTS MUST BE RESET TO TERRAIN-FOLLOWING ON 1ST
C  CALL TO AVOID NEGATIVE VALUES LEFT OVER FROM PRECEDING CALCULATIONS
C
       zrise=sfchi-sfclow
       relht=avthk
C
C  GET HEIGHT OF EACH 1ST GUESS SURFACE RELATIVE TO TERRAIN FOR EACH 
C  GRID PT. CHANGE IN HT RELATIVE TO THE LOWEST HT (MSL) IS ASSUMED 
C  PROPORTIONAL TO SIGMA FOR THAT SFC.
C
       DO 67 jy=1,nrow
          DO 66 ix=1,ncol
                 DO 65 kz=1,nlvl
            IF (kz.EQ.1) THEN
                   rhs(ix,jy,kz) = z10
            ELSE
                   rhs(ix,jy,kz) =sigma(kz)*avthk- 
     $                   (sfcht(ix,jy)-sfclow)*(1.0-slfac)
            END IF
C
                    IF (rhs(ix,jy,kz) .LE. 0.0) THEN
               WRITE (*,*) 'BAD TERRAIN-FOLLOWING IN SFCTRP'
              
               stop
            END IF
65               CONTINUE
66            CONTINUE
67     CONTINUE
C
C*********END BUG FIX WHEN MULTIPLE CASES RUN*************
C
C  INTERPOLATING WIND COMPONENTS BETWEEN SURFACE OBSERVATIONS.
C
       numtru=0
       DO 80 iob=1,numnws
          IF (ucomp(iob).GT.-9998.9) THEN
             numtru=numtru+1
             uob(numtru)=ucomp(iob)
             vob(numtru)=vcomp(iob)
             xval(numtru)=xg(iob)
             yval(numtru)=yg(iob)
          END IF    
80     CONTINUE
C
C  NOW ASSIGN INTERPOLATED VALUES TO LEVEL 1
C
       DO 124 iy=1,nrow
          y2=float(iy)
          DO 122 ix=1,ncol
             x2=float(ix)
             CALL rinvmod(trpval,x2,y2,xval,yval,numtru,uob)
         u(ix,iy,1)=trpval
             CALL rinvmod(trpval,x2,y2,xval,yval,numtru,vob)
         v(ix,iy,1)=trpval
             levbot(ix,iy)=2
C
C  ON 1ST CALL EXTRAPOLATE OR INTERPOLATE TO NEXT SURFACE
C
                  CALL lgntrp(u(ix,iy,2),z0,z10,
     $                    rhs(ix,iy,2),zero,u(ix,iy,1))
                  CALL lgntrp(v(ix,iy,2),z0,z10,
     $                    rhs(ix,iy,2),zero,v(ix,iy,1))
122       CONTINUE
C
124    CONTINUE
C
    ELSE
C
C  ON SECOND CALL USE VALUES FROM 1ST CALL TO 
C  DETERMINE COMPONENTS AT 1ST ABOVE GROUND LEVEL
C
       DO 150 iy=1,nrow
              DO 148 ix=1,ncol
C
C  NOW SET SUBSFC VALUES TO ZERO AND INTER(EXTRA)POLATE TO 
C  LOWEST ABOVE-GROUND FLOW SURFACE.
C
             IF (levbot(ix,iy).GT.2) THEN
                    DO 145 l=2,levbot(ix,iy)-1
C
                       u(ix,iy,l)=0.0
                       v(ix,iy,l)=0.0
145                 CONTINUE
         END IF
C
                 uu =u(ix,iy,1)
                 vv =v(ix,iy,1)
                 CALL lgntrp(u(ix,iy,levbot(ix,iy)),z0,z10,
     $                             rhs(ix,iy,levbot(ix,iy)),zero,uu)
                 CALL lgntrp(v(ix,iy,levbot(ix,iy)),z0,z10,
     $                             rhs(ix,iy,levbot(ix,iy)),zero,vv)
148       CONTINUE
150    CONTINUE
    END IF
C
        RETURN
C
        END
C
********************************************************************
    SUBROUTINE vrsmoo(VAL,VALTRP,NUMBER)
C*******************************************************************
C  
C  THIS PROGRAM INTERPOLATES TO FIND VALUES (VALTRP) AT GRID POINTS
C  FROM OBSERVED VALUES (VAL)
C
C       FLUDWIG, OCT 97
C  THIS VERSION USES SIMPLE INVERSE DISTANCE SQUARED WEIGHTING SO THAT
C  THOSE PARAMETERS LIKE ALTIMETER SETTING AND POTENTIAL TEMPERATURE
C  WHICH ARE NOT EXPECTED TO HAVE VERY SHARP HORIZONTAL GRADIENTS WILL
C  BE SMOOTHED.  PRESSURE OBSERVATIONS SEEMED SUFFICIENTLY UNRELIABLE
C  THAT SMOOTHING WAS DEEMED DESIRABLE.
C
C       FLUDWIG,  6/2002
C
C
C  XG,YG,VAL    COORDINATES (GRID UNITS) & VALUE AT OBSERVING SITES
C  VAL      OBSERVED VARIABLE VALUES
C  XVAL,YVAL    COORDINATES NORMALIZED FOR MQ INTERPOLATION
C  VALTRP   INTERPOLATED VALUES AT GRID POINTS
C  NUMBER   NUMBER OF OBSERVATIONS
C  
        include 'NGRIDS.PAR'
C
        include 'LIMITS.CMM'
        include 'STALOC.CMM'
C
    REAL VAL(nsites),VALTRP(nxgrd,nygrd)
    SAVE
C
C BEGIN LOOP S THROUGH GRID POINTS
C
    DO 370 ix=1,ncol
       x2=float(ix)
       DO 365 iy=1,nrow
          y2=float(iy)
          numtru=0
          sum=0.0
          sumwt=0.0
          DO 355 is=1,number
C
C CHECK THAT DATA ARE GOOD
C
             IF (val(is).GT.-9998.9) THEN
            weight=wndwt(x2,y2,xg(is),yg(is))
            sumwt=sumwt+weight
            sum=sum+weight*val(is)
         END IF
355       CONTINUE
          IF (sumwt .GT. 0.0 ) THEN
             valtrp(ix,iy)=sum/sumwt
          ELSE
             WRITE (*,*) ' NO GOOD PRESSURE OR TEMP OBS IN VRSMOO'
         WRITE (*,*) 'RETURN TO QUIT'
         pause
         stop
          END IF
365    CONTINUE
370 CONTINUE
C
    RETURN
C
    END
C
********************************************************************
    SUBROUTINE vrsdis(VAL,VALTRP,NUMBER)
C*******************************************************************
C  
C  THIS PROGRAM INTERPOLATES TO FIND VALUES (VALTRP) AT GRID POINTS
C  FROM OBSERVED VALUES (VAL)
C
C       FLUDWIG, OCT 97
C
C
C  MODIFIED TO USE INPUT VALUE WHEN WITHIN 0.5*SQRT(2) GRID UNITS
C  OR INVERSE DISTANCE POLYNOMIAL FIT BASED ON NEAREST OBSERVATIONS
C  OTHERWISE  --- RINVMOD (TRPVAL,X0,Y0,XX,YY,NOBS,VARBL)
C 
C
C       FLUDWIG,  7/2002
C
C
C  XG,YG,VAL    COORDINATES (GRID UNITS) & VALUE AT OBSERVING SITES
C  VALTRP   INTERPOLATED VALUES AT GRID POINTS
C  NUMBER   NUMBER OF OBSERVATIONS
C  
        include 'NGRIDS.PAR'
C
        include 'LIMITS.CMM'
        include 'STALOC.CMM'
C
    REAL VAL(nsites),VALTRP(nxgrd,nygrd)
    REAL VARBL(nsites),XVAL(nsites),YVAL(nsites)
C
    SAVE
C
C CHECK THAT DATA ARE GOOD
C   
    WRITE (*,*) number,(nint(100.0*val(iob)),iob=1,number)
C
    IF (number.LE.0) THEN
       WRITE(*,*) 'BAD OBSERVED DATA'
       pause
       stop
    ELSE
C
C BEGIN LOOP S THROUGH GRID POINTS
C
       DO 370 iy=1,nrow
          y2=float(iy)
          DO 365 ix=1,ncol
             x2=float(ix)
             CALL rinvmod(trpval,x2,y2,xval,yval,numtru,varbl)
         valtrp(ix,iy)=trpval
365       CONTINUE
          WRITE (*,*) iy,(nint(100.0*valtrp(jx,iy)),jx=20,40)
370    CONTINUE
    END IF
C
    RETURN
C
    END
C**********************************************************************
        REAL FUNCTION sloper(HERE,SFCMIN,HIRISE)
C**********************************************************************
C
C  DETERMINES RATIO OF HEIGHT ABOVE LOWEST TOPOGRAPHY TO MAXIMUM
C  HEIGHT ABOVE LOWEST POINT. OTHER RELATIONSHIPS CAN BE SUBSTITUTED
C  TO GIVE DIFFERENT FLOW SURFACE HEIGHTS.  -- LUDWIG 11/87
C
        IF ((hirise-sfcmin) .EQ. 0.0 .OR. hirise .EQ. 0.0) THEN
C
C  FLAT TERRAIN ALWAYS GIVES FLAT SFCS
C
           sloper =0.0
        ELSE
C
C  SFC RISE AT THIS PT PROPORTIONAL TO TERRAIN INCREMENT AS A FRACTION 
C  OF MAXIMUM TERRAIN ELEVATION DIFFERENCES.
C
           sloper = (here-sfcmin)/hirise
        END IF
C
        RETURN
        END