windest.f

Go to the documentation of this file.

C*******************************************************************
        SUBROUTINE fixwnd (IFXPT,LVL)
C******************************************************************
C
C  THIS ROUTINE IDENTIFIES PTS NEAR OBSERVATION SITES SO THAT
C  WIND ADJUSTMENTS CAN BE RESTRAINED IN SUBROUTINE BAL5.
C
C___________________________________________________________________
C
C  MODIFIED SO THAT ALL LEVELS ABOVE AN UPPER WIND SITE ARE IDENTIFIED,
C  BUT ONLY THE FIRST ABOVE GROUND LEVEL FOR SURFACE WIND SITES.
C
C       F. L. LUDWIG, 3/96
C___________________________________________________________________
C
C  FURTHER MODIFIED SO THAT ABOVE GROUND POINTS ON A SURFACE ARE 
C  FLAGGED FOR RESTRAINED ADJUSTMENT IF THEY HAVE 3 OR MORE OF THE 4 
C  SURROUNDING POINTS ARE BELOW GROUND.
C
C       FLUDWIG 3/2000
C___________________________________________________________________
C
          include 'NGRIDS.PAR'       
C
          include 'FLOWER.CMM'
          include 'LIMITS.CMM'
          include 'STALOC.CMM'
C
        LOGICAL IFXPT(nxgrd,nygrd)
    SAVE
C
C  INTIALIZE IFXPT VALUES TO FALSE
C
    DO 22 ix=1,nxgrd
       DO 20 iy=1,nygrd
        ifxpt(ix,iy)=.false.
20     CONTINUE
22  CONTINUE
C
C  CHECK ALL GRID POINTS TO SEE IF THEY SHOULD BE ADJUSTED
C  IF THE POINT ITSELF, OR 3 OR 4 OF THE SURROUNDING POINTS 
C  ARE SUBTERRAINIAN, RESTRICT ADJUSTMENTS, I.E. SET IFXPT TO .TRUE.  
C
    DO 50 ix=2,ncol-1
       DO 48 iy=2,nrow-1
          IF (rhs(ix,iy,lvl).LE. 0.0) THEN
         ifxpt(ix,iy)=.true.
          ELSE
         nsubter=0
         DO  44 jx=ix-1,ix+1,2
            DO  40 jy=iy-1,iy+1,2
                IF (rhs(jx,jy,lvl).LE. 0.0) 
     $                             nsubter=nsubter+1
40          CONTINUE
44      CONTINUE
        IF (nsubter .GE. 3)  THEN
                   ifxpt(ix,iy)=.true.
        ELSE
                   ifxpt(ix,iy)=.false.
        END IF
         END IF
48    CONTINUE
50  CONTINUE
              
C
C  SET FLAG FOR LIMITED ADJUSTMENT AT POINTS AROUND OBS SITE
C
    DO 200 i=1,nsites
          IF (i .LE. numnws ) THEN
C
C  CHECK TO SEE IF THIS SITE HAD AN OBSERVATION
C
          IF (jgood(i)) THEN
                  ix=nint(xg(i))
                  iy=nint(yg(i))
                  IF (ix.GT.0 .AND. ix.LE.ncol   
     $                 .AND. iy.GT.0 .AND.    
     $                          iy .LE.nrow) THEN
C
C  CHECK TO SEE IF THIS IS FIRST LEVEL ABOVE SURFACE FOR 
C  SURFACE OBSERVATIONS.  
C
             IF (rhs(ix,iy,lvl) .GT. 0.0 .OR. lvl.EQ.1) 
     $                               ifxpt(ix,iy)=.true.
          END IF
           END IF
        END IF
C
200 CONTINUE
C
        RETURN
        END
C
C*********************************************************************
        SUBROUTINE betwin
C*********************************************************************
C
C  THIS SUBROUTINE ESTIMATES WINDS BETWEEN THE TOP AND BOTTOM LEVELS. 
C  THE DEVIATION OF OBSERVED WINDS FROM A LOG PROFILE IS FIRST 
C  DETERMINED.  THEN THE DEVIATIONS AT EACH LEVEL ARE INTERPOLATED BY 
C  AN INVERSE DISTANCE TO A POWER WEIGHTING SCHEME.  THE 
C  INTERPOLATED DEVIATIONS ARE USED TO CORRECT  THE CALCULATED LOG 
C  PROFILES AT THE GRID POINTS.
C               --F  LUDWIG  12/87
C
    include 'NGRIDS.PAR'
C
    include 'ANCHOR.CMM'
    include 'FLOWER.CMM'
    include 'LIMITS.CMM'
    include 'STALOC.CMM'
    include 'TSONDS.CMM'
C
        REAL DDOPU(nwsite,nzgrd),DDOPV(nwsite,nzgrd)
        REAL DUTMP(nsites),DVTMP(nsites),XVAR(nsites),YVAR(nsites)
C
    INTEGER LBOT(nwsite)
C
    LOGICAL OKSOND(nwsite,nzgrd)
C
        IF (nlvl .LE. 3) THEN
           WRITE (*,*) ' ONLY ',nlvl-1,' FLOW SURFACES'
           RETURN
        END IF
C
C  GET LOG PROFILES AT OBSERVATION POINTS AND DEVIATIONS FROM THEM. 
C  Z0=ROUGHNESS HT.
C
    IF (numdop .GT. 0) THEN
       DO 50 jdop = 1,numdop
C
C  FIND LOWEST FLOW SURFACE FOR THIS SONDE
C
          DO 22 il=2,nlvl
             IF (zsigl(jdop,il) .GE. z0) THEN
                    h0=zsigl(jdop,il)
            lbot(jdop)=il
            GO TO 24
         END IF
22        CONTINUE
C
C  INTERPOLATE BETWEEN LOWEST AND HIGHEST LEVELS ABOVE THE
C  SURFACE FOR GETTING DEVIATIONS.
C
24        uu = usig(jdop,lbot(jdop))
          vv = vsig(jdop,lbot(jdop))
          h0 = zsigl(jdop,lbot(jdop))
              utop = usig(jdop,nlvl)
              vtop = vsig(jdop,nlvl)
              ztop = zsigl(jdop,nlvl)
              ddopu(jdop,nlvl)=0.0
              ddopv(jdop,nlvl)=0.0
          oksond(jdop,nlvl)=.true.
          oksond(jdop,lbot(jdop))=.true.
              DO 38 ll = 2,nlvl-1
                 IF (ll .GT. lbot(jdop)) THEN
                    zz = zsigl(jdop,ll)
                    CALL lgntrp(uu,h0,ztop,zz,u0,utop)
                    CALL lgntrp(vv,h0,ztop,zz,v0,vtop)
C
C  GET DEVIATIONS FROM LOG INTERPOLATED VALUE FOR THIS LEVEL
C
                    ddopu(jdop,ll) = usig(jdop,ll)-uu
                    ddopv(jdop,ll) = vsig(jdop,ll)-vv
                oksond(jdop,ll)=.true.
                 ELSE
                    ddopu(jdop,ll)=0.0
                    ddopv(jdop,ll)=0.0
                IF (ll.NE.lbot(jdop)) oksond(jdop,ll)=.false.
                 END IF
38            CONTINUE
50         CONTINUE    
    END IF
C
C  GET LOG PROFILES AT EACH GRID POINT-- FROM 1ST ABOVE-GROUND LEVEL.
C  THEN ADD DEVIATION FROM LOG AS DETERMINED FROM HORIZONTAL 
C  INTERPOLATION BETWEEN SOUNDINGS
C
    nevent=0
        DO 100 ix = 1,ncol
       xhere=float(ix)
           DO 90 iy = 1,nrow
          yhere=float(iy)
              utop = u(ix,iy,nlvl)
              vtop = v(ix,iy,nlvl)
              ztop = rhs(ix,iy,nlvl)
              h0 = 10.0
              u0 = u(ix,iy,1)
              v0 = v(ix,iy,1)
          DO 80 ll=levbot(ix,iy),nlvl-1
C
C  GET COMPONENTS FROM LOG LINEAR INTERPOLATION BETWEEN TOP AND BOTTOM
C
                 zz=rhs(ix,iy,ll)
                 CALL lgntrp(u(ix,iy,ll),h0,ztop,zz,u0,utop)
                 CALL lgntrp(v(ix,iy,ll),h0,ztop,zz,v0,vtop)
C
C  GET CORRECTION FOR THIS GRID POINT AND FLOW LEVEL
C
             n4use=0
             DO 78 jdop = 1,numdop
C
            IF (oksond(jdop,ll)) THEN
                   n4use=n4use+1
               dutmp(n4use)=ddopu(jdop,ll)
               dvtmp(n4use)=ddopv(jdop,ll)
               xvar(n4use)=xdop(jdop)
               yvar(n4use)=ydop(jdop)
             END IF
78            CONTINUE
C
C  IF THIS GRID PT. IS BELOW THE BOTTOM FLOW SFC FOR ALL SOUNDINGS, WE 
C  SET THE DEVIATION TO 0 AND USE THE INTERPOLATED VALUE
C
              IF (n4use.LE.0) THEN
             nevent=nevent+1
             du=0.0
                 dv=0.0
          ELSE
                 CALL rinvmod(du,xhere,yhere,xvar,yvar,
     $                             n4use,dutmp)
                 CALL rinvmod(dv,xhere,yhere,xvar,yvar,
     $                             n4use,dvtmp)
          END IF
C
C  CHECK FOR INTERPOLATED SOND DATA AND USE IT TO CORRECT  THE
C  LOG PROFILE.I IF NO OBS ON THIS FLOW SFC, USE THE VALUES ALREADY
C  INTERPOLATED.
C
                  u(ix,iy,ll)=du+u(ix,iy,ll)
                  v(ix,iy,ll)=dv+v(ix,iy,ll)
C              
80             CONTINUE
90     CONTINUE
100 CONTINUE
C
C  INTRODUCE INFLUECE OF OBSERVED SURFACE WINDS ON LOWEST LAYER ALOFT 
C  BY INTERPOLATING BETWEEN ANEMOMETER HT. AND 2ND LOWEST SFC.TO 
C  OBTAIN WIND ON THE LOWEST ABOVE GROUND SURFACE.   IF ONLY THE TOP
C  SFC IS ABOVE GROUND, EXTRAPOLATE UP FROM SFC.   
C
    DO 200 ix=1,ncol
       DO 190 iy=1,nrow
C
          IF (levbot(ix,iy).LT.nlvl) THEN
             u0=u(ix,iy,1)
             v0=v(ix,iy,1)
             h0=z10
         utop=u(ix,iy,levbot(ix,iy)+1)
         vtop=v(ix,iy,levbot(ix,iy)+1)
         ztop=rhs(ix,iy,levbot(ix,iy)+1)
         zz=rhs(ix,iy,levbot(ix,iy))
                 CALL lgntrp(uu,h0,ztop,zz,u0,utop)
                 CALL lgntrp(vv,h0,ztop,zz,v0,vtop)
         u(ix,iy,levbot(ix,iy))=uu
         v(ix,iy,levbot(ix,iy))=vv
          ELSE
         u0=0.0
         v0=0.0
         h0=z0
         utop=u(ix,iy,1)
         vtop=v(ix,iy,1)
         ztop=z10
         zz=rhs(ix,iy,nlvl)
                 CALL lgntrp(uu,h0,ztop,zz,u0,utop)
                 CALL lgntrp(vv,h0,ztop,zz,v0,vtop)
         u(ix,iy,levbot(ix,iy))=uu
         v(ix,iy,levbot(ix,iy))=vv
          END IF
190    CONTINUE
200 CONTINUE
C
        RETURN
C
        END