richardson.f

Go to the documentation of this file.

C******************************************************************
    SUBROUTINE getrich
C******************************************************************
C
C  GETS NECESSARY WIND SHEAR AND POTENTIAL TEMPERATURE GRADIENTS AT
C  EACH GRID POINT FROM DIFFERENCES AND INTERPOLATION ON EACH SURFACE.
C  FROM THESE THE RICHARDSON NUMBER IS ESTIMATED.  VALUES OF POTENT.
C  TEMP. AND PRESSURE ARE ALSO ESTIMATED
C
C   FLUDWIG 3/02
C
    include 'NGRIDS.PAR'
C
    include 'ANCHOR.CMM'
    include 'FLOWER.CMM'
    include 'LIMITS.CMM'
        include 'STALOC.CMM'
    include 'TSONDS.CMM'
C
    REAL HTDIF,TMPDIF,TMPAVG,SQDIF
C
C  WE NOW HAVE TOP AND BATTOM EVERYWHERE VALUES FOR ALTIMETER SETTING
C  AND POTENTIAL TEMPERATURE.  NOW FILL IN BETWEEN VALUES.
C
    DO 180 ix=1,nxgrd
       DO 160 iy=1,nygrd
               DO 140 iz=1,nflat
C
C  IDENTIFY GRID POINTS AT OR BELOW 20 M ABOVE SFC AND FLAG HEIGHT
C  DIFFERENCE AS NEGATIVE FOR FUTURE SKIPPING.  GET DIFFERENCE IN 
C  HEIGHT BETWEEN SURFACES TO BE USED FOR CALCULATING BRUNT VAISALA 
C  AND RICHARDSON VALUES
C
        zblow=100.0
             zatxy=sfcht(ix,iy)+2.0*z10
             IF (iz.EQ.1)  THEN
            htdif=-9.0
         ELSE IF (iz.EQ.2) THEN
            zblow=zchooz(2)
            htdif=zchooz(3)-zchooz(2)
         ELSE IF (iz.GT.2 .AND. iz.LT.nflat) THEN
            zblow=zchooz(iz-1)
            htdif=zchooz(iz+1)-zchooz(iz-1)
         ELSE IF (iz.EQ.nflat) THEN
            zblow=zchooz(nflat-1)
            htdif=zchooz(nflat)-zchooz(nflat-1)
         END IF
C
C  FLAG SFC BELOW GROUND
C
         IF (doflat .AND. zblow .LT. zatxy) htdif=-9.
C
C  GRADIENTS UNRELIABLE AT ANEMOMETER LEVEL.  SET BULK RICHARDSON 
C  NUMBER & BRUNT VAISALA PERIOD AS MISSING WHEN FLAGGED HTDIF < 0.
C
             IF (htdif .LT. 2.0*z10) THEN
                irngrf(ix,iy,iz)=-2
                ibvgrf(ix,iy,iz)=0
         ELSE
            IF (iz.EQ.nflat) THEN
               tmpdif=0.1*float(iptgrf(ix,iy,nflat)-
     $                                  iptgrf(ix,iy,nflat-1))
               tmpavg=0.05*float(iptgrf(ix,iy,nflat)+
     $                                  iptgrf(ix,iy,nflat-1))
               iudif=iugraf(ix,iy,nflat)-iugraf(ix,iy,nflat-1)
               ivdif=ivgraf(ix,iy,nflat)-ivgraf(ix,iy,nflat-1)
            ELSE IF (iz.EQ.2) THEN
               tmpdif=0.1*float(iptgrf(ix,iy,3)-
     $                                         iptgrf(ix,iy,2))
               tmpavg=0.05*float(iptgrf(ix,iy,3)+
     $                                         iptgrf(ix,iy,2))
               iudif=iugraf(ix,iy,3)-iugraf(ix,iy,2)
               ivdif=ivgraf(ix,iy,3)-ivgraf(ix,iy,2)

            ELSE 
               tmpdif=0.1*float(iptgrf(ix,iy,iz+1)-
     $                                      iptgrf(ix,iy,iz-1))
               tmpavg=0.05*float(iptgrf(ix,iy,iz+1)+
     $                                      iptgrf(ix,iy,iz-1))
               iudif=iugraf(ix,iy,iz+1)-iugraf(ix,iy,iz-1)
               ivdif=ivgraf(ix,iy,iz+1)-ivgraf(ix,iy,iz-1)
            END IF
C
C  SQUARE SPEED DIFFERENCE & CONVERT FROM (CM/S)**2 TO (M/S)**2
C
                    sqdif=float(iudif**2)+float(ivdif**2)
            sqdif=1.0e-4*sqdif
C
C  NOW GET BRUNT VAISALA PERIOD AND BULK RICHARDSON NUMBER
C
               ibvgrf(ix,iy,iz)=nint(bvperd(htdif,tmpdif,tmpavg))
               IF (ibvgrf(ix,iy,iz) .GT. 900)  
     $                              ibvgrf(ix,iy,iz)=900
               irngrf(ix,iy,iz)=
     $                   nint(1000.0*bulkri(htdif,tmpdif,tmpavg,sqdif))
         END IF
C
140        CONTINUE
160    CONTINUE
180 CONTINUE
C
    RETURN
C
    END
C
C******************************************************************
    REAL FUNCTION bvperd (HTDIF,TMPDIF,TMPAVG)
C******************************************************************
C
C  RETURNS INVERSE OF BRUNT-VASALA FREQUENCY, CALCULATED FROM
C
C   TMPAVG  = AVERAGE POTENT. TEMP IN LAYER
C   TMPDIF  = POTENTIAL TEMPERATURE DIFFERENCE OVER HTDIF
C   HTDIF   = THICKNESS OF LAYER
C   GRAV    = GRAVITATIONAL ACCELERATION
C
C  POTENTIAL TEMP DIFFERENCE MUST BE > 0.  IF LESS, SET PERIOD
C  TO 999 SEC.
C
C   FLUDWIG 10/02
C
    parameter(grav=9.8)
C
    REAL HTDIF,TMPDIF,TMPAVG,BV2
C
    IF (tmpdif.GT. 0.0) THEN
       bv2=htdif*tmpavg/(grav*tmpdif)
       bvperd=sqrt(bv2)
    ELSE
       bvperd=999.
    END IF
C
    RETURN
C
    END
C
C******************************************************************
    REAL FUNCTION bulkri(HTDIF,TMPDIF,TMPAVG,SQDIF)
C******************************************************************
C
C  RETURNS BULK RICHARDSON NUMBER, CALCULATED FROM
C
C   TMPAVG  = AVERAGE POTENT. TEMP IN LAYER
C   SQDIF   = DU**2 + DV**2 -- SQUAREDWND DIFFERENCE OVER HTDIF
C   TMPDIF  = POTENTIAL TEMPERATURE DIFFERENCE OVER HTDIF
C   HTDIF   = THICKNESS OF LAYER
C   GRAV    = GRAVITATIONAL ACCELERATION
C
C   -0.99 RETURNED FOR ZERO WIND SHEAR.
C
C   FLUDWIG 10/02
C
    parameter(grav=9.8)
C
    REAL HTDIF,TMPDIF,TMPAVG,SQDIF
C
    IF (tmpavg*sqdif .GT. 0.0) THEN
       bulkri=(grav*tmpdif*htdif)/(tmpavg*sqdif)
       IF (bulkri .GT. 99.0) bulkri=99.0
    ELSE
       bulkri=-0.002
    END IF
C
    RETURN
C
    END