levwnd.f

Go to the documentation of this file.

C*********************************************************************
        SUBROUTINE levwnd
C*********************************************************************
C
C  INTERPOLATES MASS ADJUSTED FIELD TO ANEMOMETER HEIGHT 
C  (Z10, FOR Z INDEX=1)
C  & TO NHORIZ-1 FLAT PLANES (FOR Z INDICES 2-NHORIZ) SET AT THE 
C  HEIGHTS ZCHOOZ ABOVE THE LOWEST TERRRAIN GRID POINT IN THE COARSE 
C  GRID, AS DEFINEDTHE VARIABLE ARRAY SIGMA.  THIS SUBROUTINE ALSO 
C  CONVERTS INTERPOLATED WINDS BACK TO METEOROLOGICAL SPEEDS
C  AND ANGLES THAT CAN BE PLOTTED IF DESIRED.  LUDWIG--JANUARY 1988
C
C  MODIFIED APRIL 2002 TO INTERPOLATE OTHER VARIABLES (RICHARDSON NO.,
C  POT. TEMP. LAPSE RATE, PRESSURE, POT. TEMP, BRUNT VAISALA PERIOD) 
C  TO FLAT SURFACES AS WELL.  ALSO SIMPLIFIED LOGIC.
C
C   FLUDWIG 4/02
C
        parameter(rad2d=180./3.14159)
C
        include 'NGRIDS.PAR'
C
        include 'ANCHOR.CMM'
        include 'FLOWER.CMM'
        include 'LIMITS.CMM'
        include 'TSONDS.CMM'
C
        LOGICAL IFXPT(nxgrd,nygrd),DOTRPO,DIDFND
C
C  MAKE SURE EVERYTHING IS STILL HERE ON SUBSEQUENT CALLS
C
        SAVE
C
    k1=1
    z0=zzero
    z10=10.0
    nbada=0
C
        CALL fixwnd (ifxpt,1)
C
    DO 400 ix=1,ncol
           DO 390 iy=1,nrow
C
C  1ST LEVEL REPRESENTS A TERRAIN FOLLOWING SURFACE, Z10 M ABOVE SFC
C  EVERYWHERE.  IN ORDER TO REFLECT BOTH THE OBSERVATIONS USED TO
C  GET FIRST ESTIMATES OF ANEMOMETER LEVEL WINDS AND THE ADJUSTED 
C  WINDS ON THE FIRST ABOVE GROUND FLOW SURFACE, WE USE FOR THE LOWER
C  END OF THE INTERPOLATION, HALF THE INITIAL LEVEL 1 WIND ESTIMATE AT
C  THE HT WHERE A LOG PROFILE WOULD HAVE THE WINDS BE 50% OF THEIR 
C  MAGNITUDE AT Z=Z10.  VALUES AT THE 1ST LEVEL ABOVE GROUND ARE USED
C  FOR THE UPPER END OF THE INTERPOLATION.  FOR PTS NEAR OBS, THE 
C  WEIGHTING OF ANEMOMETER HT VALUES IS INCREASED BY USING Z FOR
C  90% OF ANEMOMETER HT VLAUE.  LAPSE RATE UNRELIABLE AT THIS LEVEL, SO
C  BULK RICHARDSON NUMBER AND BRUNT-VAISALA PERIOD NOT CALCULATED
C
          zagl=z10
          IF (ifxpt(ix,iy)) THEN
                 bottom=z0*((z10/z0)**0.9)
                 u0=0.9*u(ix,iy,1)
                 v0=0.9*v(ix,iy,1)
                 w0=0.9*w(ix,iy,1)
          ELSE
                 bottom=z0*((z10/z0)**0.5)
                 u0=0.5*u(ix,iy,1)
                 v0=0.5*v(ix,iy,1)
                 w0=0.5*w(ix,iy,1)
          END IF
C
C  NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION.
C
          top=rhs(ix,iy,levbot(ix,iy))
              u1=u(ix,iy,levbot(ix,iy))
              v1=v(ix,iy,levbot(ix,iy))
              w1=w(ix,iy,levbot(ix,iy))
C
C  WE NOW HAVE VALUES ABOVE AND BELOW THE ANEMOMETER HEIGHT -- 
C  INTERPOLATE
C
              CALL lgntrp(uuuf,bottom,top,zagl,u0,u1)
          u(ix,iy,1)=uuuf
              CALL lgntrp(vvvf,bottom,top,zagl,v0,v1)
          v(ix,iy,1)=vvvf
              CALL lgntrp(wwwf,bottom,top,zagl,w0,w1)
          w(ix,iy,1)=wwwf
              iugraf(ix,iy,1)=nint(100.0*u(ix,iy,1))
              ivgraf(ix,iy,1)=nint(100.0*v(ix,iy,1))
              iwgraf(ix,iy,1)=nint(100.0*w(ix,iy,1))
          iptgrf(ix,iy,1)=nint(10.0*tmpkel(ix,iy,1))
C         IPRGRF(IX,IY,1)=NINT(100.0*SETALT(IX,IY,1))
          patsfc=cvt2p(setalt(ix,iy,1),sfcht(ix,iy)+z10)
          iprgrf(ix,iy,1)=nint(10.0*patsfc)
C
C  DONE WITH ANEMOMETER HEIGHT VALUES, NOW INTERPOLATE TO HORIZONTAL 
C  PLANES (IZ>.
C 
          DO 300 iz=2,nflat
C  
C  GET HEIGHT OF HORIZONTAL PLANE ABOVE LOCAL GROUND LEVEL (ZAGL), THEN
C  FIND 1ST FLOW SURFACE ABOVE THIS FLAT SURFACE, MEANWHILE FILLING
C  BELOW SFC LEVELS WITH MISSING DATA VALUES.  IF ALL SURFACES ARE 
C  LOWER THAN TERRAIN AT THIS GRID POINT SET ALL LEVELS TO MISSING 
C  VALUES AND PROCEED TO NEXT GRID PT.
C
                 IF (doflat) THEN
            zagl=zchooz(iz)-sfcht(ix,iy)
            zmsl=zchooz(iz)
         ELSE
            zagl=zchooz(iz)
            zmsl=zchooz(iz)+sfcht(ix,iy)
         END IF
C
             IF (zagl.LE.z0) THEN
            dotrpo=.false.
C
C  THIS SURFACE BELOW GROUND;  SET VALUES TO ZERO OR MISSING
C       
                    iugraf(ix,iy,iz)=0
                    ivgraf(ix,iy,iz)=0
                    iwgraf(ix,iy,iz)=0
                iptgrf(ix,iy,iz)=0
                iprgrf(ix,iy,iz)=0
C
             ELSE IF (zagl .GE. rhs(ix,iy,nlvl)) THEN
            dotrpo=.false.
C
C  EXRAPOLATION SURFACE ABOVE HIGHEST FLOW SURFACE.  USE VALUES
C  BASED ON THOSE ON HIGHEST SURFACE. ASSUME A LAPSE RATE OF 
C
                    iugraf(ix,iy,iz)=nint(100.0*u(ix,iy,nlvl))
                    ivgraf(ix,iy,iz)=nint(100.0*v(ix,iy,nlvl))
                    iwgraf(ix,iy,iz)=nint(100.0*w(ix,iy,nlvl))
            IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
               ttt=tmpkel(ix,iy,nlvl)+dthdzl(ix,iy,nlvl)*
     $                                     (zagl-rhs(ix,iy,nlvl))
               iptgrf(ix,iy,iz)=nint(10.0*ttt)
               atop=setalt(ix,iy,nlvl)
C              IPRGRF(IX,IY,IZ)=NINT(100.0*ATOP)
               iprgrf(ix,iy,iz)=nint(10.0*cvt2p(atop,zmsl))
            END IF
C
         ELSE IF (zagl.LE. z10) THEN
            dotrpo=.true.
C
C  INTERPOLATION SURFACE IS BELOW ANEMOMETER HEIGHT.  USE LOG PROFILE 
C  WIND VALUES FROM ANEMOMETER LEVEL.  USE ANEMOMETER LEVEL VALUES FOR 
C  OTHER PARAMETERS.
C 
            ibr=4
            bottom=z0
            u0=0.0
            v0=0.0
            w0=0.0
            IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
               t0=tmpkel(ix,iy,1)
               t1=t0
C
C ASSUME AT LOWEST SURFACE  APPLIES HERE.
C
               alt0=setalt(ix,iy,levbot(ix,iy))
C              IPRGRF(IX,IY,IZ)=NINT(100.0*ALT0)
               iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
            END IF
C
C  NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION.  ASSUME NO 
C  CHANGE AND USE SAME VALUES FOR THE VARIABLES OTHER THAN WIND.
C 
            top=z10
            u1=u(ix,iy,1)
            v1=v(ix,iy,1)
            w1=w(ix,iy,1)
C
             ELSE IF (zagl .LT.  rhs(ix,iy,levbot(ix,iy))
     $                                     .AND. zagl.GT.z10) THEN
            dotrpo=.true.
C
C  THIS INTERPOLATION SFC MUST BE BETWEEN ANEMOMETER AND THE
C  LOWEST ABOVE-GROUND FLOW SURFACE.
C
            ibr=2
C
C  GET VALUES TO USE FOR BOTTOM END OF INTERPOLATION.
C
            bottom=z10
                    u0=u(ix,iy,1)
                    v0=v(ix,iy,1)
                    w0=w(ix,iy,1)
            IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
               t0=tmpkel(ix,iy,1)
               t1=tmpkel(ix,iy,levbot(ix,iy))
               alt0=setalt(ix,iy,levbot(ix,iy))
C              IPRGRF(IX,IY,IZ)=-NINT(100.0*ALT0)
               iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
            END IF
C
C  NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION.
C
                    top=rhs(ix,iy,levbot(ix,iy))
                    u1=u(ix,iy,levbot(ix,iy))
                    v1=v(ix,iy,levbot(ix,iy))
                    w1=w(ix,iy,levbot(ix,iy))
C
         ELSE IF (zagl.GE. rhs(ix,iy,levbot(ix,iy)) 
     $                        .AND. zagl .LT. rhs(ix,iy,nlvl)) THEN
            dotrpo=.true.
C
C  GETTING VALUES WHEN INTERPOLATION SFC IS BETWEEN THE LOWEST AND 
C  THE HIGHEST FLOW SFCS
C
            DO 130 ll=levbot(ix,iy),nlvl-1
               didfnd=.false.
               IF (zagl.GE.rhs(ix,iy,ll) .AND.
     $                                 zagl.LT.rhs(ix,iy,ll+1)) THEN                                
C
C  FOR INTERPOLATION SFC BETWEEN 2 FLOW SFCS
C
                  ibr=3
                          bottom=rhs(ix,iy,ll)
                          u0=u(ix,iy,ll)
                          v0=v(ix,iy,ll)
                          w0=w(ix,iy,ll)
C
C  WIND VALUES TO USE FOR TOP END OF INTERPOLATION.
C
                          top=rhs(ix,iy,ll+1)
                          u1=u(ix,iy,ll+1)
                          v1=v(ix,iy,ll+1)
                          w1=w(ix,iy,ll+1)
C
C  USE WEIGHTED ALTIMETER SETTING FROM FLOW SFCS.
C
                      IF (zagl.EQ.rhs(ix,iy,ll)) THEN
                             alt0=setalt(ix,iy,ll)
              ELSE IF (rhs(ix,iy,ll+1) .EQ. zagl) THEN
                 alt0=setalt(ix,iy,ll+1)
              ELSE
                         wt1=1.0/abs(zagl-rhs(ix,iy,ll))
                         wt2=1.0/abs(zagl-rhs(ix,iy,ll+1))
                 alt0=(wt1*setalt(ix,iy,ll)+
     $                              wt2*setalt(ix,iy,ll+1))/(wt1+wt2)
              END IF
C                 IPRGRF(IX,IY,IZ)=NINT(100.0*ALT0)
                  iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
              IF (dothet .OR. dopres .OR. dobri
     $                                       .OR. dobvpd) THEN
                 t0=tmpkel(ix,iy,ll)
                 t1=tmpkel(ix,iy,ll+1)
              END IF
C
C WE FOUND WHERE IT IS, WE CAN GO ON TO THE NEXT STEP
C
              didfnd=.true.
                   END IF
               IF (didfnd) GO TO 135
130             CONTINUE
135             CONTINUE
C
         END IF
C
C  IF WE HAVE NOT ALREADY DEFINED THE VALUES, THEN WE NEED 
C  TO INTERPLATE, USING VALUES DEFINED ABOVE.
C
         IF (dotrpo) THEN
C
C  IF ONE BVP VALUE, GOOD USE IT
C
                    CALL lgntrp(uuuf,bottom,top,zagl,u0,u1)
                    CALL lgntrp(vvvf,bottom,top,zagl,v0,v1)
                    CALL lgntrp(wwwf,bottom,top,zagl,w0,w1)
                iugraf(ix,iy,iz)=nint(100.0*uuuf)
                ivgraf(ix,iy,iz)=nint(100.0*vvvf)
                iwgraf(ix,iy,iz)=nint(100.0*wwwf)
C
C  VARIABLES ARE SCALED TO INTEGERS FOR WRITING FILES.  UNITS ARE
C  CM/SEC FOR WIND COMPONENTS, (DEG K)/10 FOR POTENTIAL TEMP., MB (HP) 
C  FOR PRESSURE, SECONDS FOR BRUNT VAISALA PERIOD, AND 100* BULK 
C  RICHARDSON NO.
C
            aaset=xlintr(alt0,alt1,zagl,bottom,top)
            IF (dothet .OR. dopres .OR. dobri.OR. dobvpd) 
     $                  iptgrf(ix,iy,iz)=nint(10.0*
     $                              xlintr(t0,t1,zagl,bottom,top))          

         END IF
C
300       CONTINUE
C
390    CONTINUE
400 CONTINUE
C
C  CONVERT COMPONENTS IN HORIZONTAL OR TERRAIN FOLLOWING 
C  PLANES TO METEOROLOGICAL WINDS SPEED & DIRECTION 
C
    DO 500 lev=1,nflat
           DO 490 ix=1,ncol
          DO 480 iy=1,nrow
C
C   CONVERT INTEGER CM/S BACK TO REAL M/S
C
                 ug=0.01*float(iugraf(ix,iy,lev))
                 vg=0.01*float(ivgraf(ix,iy,lev))
C
                 IF (ug.EQ.zero.AND.vg.EQ.zero) THEN
                    spdmet(ix,iy,lev)=zero
                    dirmet(ix,iy,lev)=zero
                 ELSE
                    spdmet(ix,iy,lev)=sp(ug,vg)
                    dirmet(ix,iy,lev)=dd(ug,vg)
                 END IF
480       CONTINUE
490    CONTINUE
500 CONTINUE       
C
    RETURN     
C
    END