      SUBROUTINE WVARP(rhobar,omega2,radei,drag,hcorfe,
     &                IREG,NWF,SINY,COSY,WLAT,WLON,WSTD,WDRG,
     &                BC,BF,SCALE)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  NASA / DAO Optimum Interpolation Assimilation System - Version 1.3
C
C  Program Name: WVARP
C  Date        : May, 1992
C  Programmer  : S. Bloom, NASA/SDUO
C
C WVAR FORMS THE ARRAY OF WIND FCST ERROR VARIANCES (IN CORRECT UNITS
C
C  INPUT ARGUMENTS:
C
C     NWF      = NW + NH
C     WSTD     = PRESSURE VARIANCES AT GRID POINTS
C
C  OUTPUT ARGUMENTS:
C
C                       RHO * F
C     BC       = ----------------------  ("CORIOLIS TYPE" SFC COEF)
C                (RHO*F)**2 + (DRAG)**2
C
C                        DRAG
C     BF       = ----------------------  ("DRAG TYPE" SFC COEF)
C                (RHO*F)**2 + (DRAG)**2
C
C     SCALE    = DIMENSIONAL FACTOR FROM CORRELATION TABLE
C     WSTD     = WIND ERROR VARIANCE (IN UNITS OF M SEC-1)
C
C
C
C ===============================================================
C
C     TROPICAL DECOUPLING BUILT INTO WSTD USING RELATION
C
C
C             WSTD = (1-A) * TWSTD  +  (A) * WSTD(EKMAN)
C
C  --------------------------------------------------------------
C
C                                   0       ,    ABS(LAT) <= 10
C
C   WITH 'A' A LINEAR RAMP:   .1*ABS(LAT)-1 , 10<ABS(LAT) <= 20
C
C                                   1       ,    ABS(LAT) >  20
C
C
C
C     FOR CURRENT TESTS...USING TWSTD = 5 M SEC-1
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      real rhobar,omega2,radei
      real hcorfe(1000,3)
      DIMENSION BC(700),BF(700),COSY(700),SINY(700)
      dimension WLAT(700),WLON(700),WSTD(700)
      DIMENSION WDRG(700)


      DATA ZER0/0.0/
      DATA ONE/1.0/
      DATA TWO/2.0/
      DATA FIVE/5.0/
      DATA TEN/10.0/
      DATA TWENTY/20.0/
      DATA TENTH/0.1/
      DATA TENM5/1.E-5/
      data zero / 0.0 /

c      CALL GETCONH('RHOBAR  ',RHOBAR)
c      CALL GETCONH('OMEGA2  ',OMEGA2)
c      CALL GETCONH('RADEI   ',RADEI)
C
C.... OBTAIN SCALE FACTOR FROM CORRELATION TABLE
C

      ICOR=1

c      SCALE = HCORFE(1,3,IREG,ICOR)
      SCALE = HCORFE(1,3)

c      IF(SCALE .GE. ZERO) then
c        message = ' WVARP: Bad scale value detected. '
c        call flash( 5,message,1 )
c      end if

      SCALE =  SQRT(-SCALE)

c      CALL TERPV(WLAT,WLON,WDRG,NWF,DRAG)

C
C..... TEMPORARY (?) CODE TO ASSIGN GLOBAL CONSTANT VALUE TO WDRG
C-----------------------------
      DO 3 I=1,NWF
              WDRG(I) = TENM5
 3    CONTINUE
C-----------------------------

      DO 5 I=1,NWF

              WDRG(I) = RHOBAR * WDRG(I)

 5    CONTINUE

C
C.... SET UP SURFACE ANALYSIS EKMAN BALANCE COEFFCIENTS
C

      DO 10 I=1,NWF

              RF    =  RHOBAR*OMEGA2*SINY(I)
              BC(I) =  TENTH/(RF**2 + WDRG(I)**2)

 10   CONTINUE

      DO 12 I=1,NWF

              BF(I) = WDRG(I) * BC(I)

 12   CONTINUE

      DO 15 I=1,NWF

              BC(I) = RHOBAR * OMEGA2 * SINY(I) * BC(I)

 15   CONTINUE

C
C.... HERE WE CONVERT WSTD INTO A EKMAN-BALANCE ESTIMATE OF WIND
C     ERROR VARIANCE...USING FORMULA:
C
C
C
C        WSTD(EKMAN) = SCALE * SQRT(BC**2 + BF**2) * WSTD(OLD) /RADE
C
C     NOTE, THE UNITS OF WSTD(NEW) ARE M SEC-1
C     AND THAT WE INITIALLY PUT EKMAN BALANCE ESTIMAE IN WSTD FOR ALL
C     LATITUDES
C

      DO 20 I=1,NWF

              BFAC    = BC(I)**2 + BF(I)**2
              BFACQ   = SQRT(BFAC)

              WSTD(I) = SCALE * BFACQ * WSTD(I) * RADEI

              BC(I)   = BC(I)/BFACQ
              BF(I)   = BF(I)/BFACQ


 20   CONTINUE

C
C.... NOW, WE WANT TO PUT TROPICAL DECOUPLING INTO WSTD:
C
C
C        = TWSTD                       , FOR ABS(LAT) LE 10
C
C WSTD   = (1-A)*TWSTD + (A)*WSTD(EK)  , FOR 0 LT ABS(LAT-10) LE 10
C
C        = WSTD(EK)                    , FOR ABS(LAT) GT 20
C
C
C     HERE 'A' IS A SIMPLE LINEAR RAMP FUNCTION GOING FROM 0 AT 10
C     DEGREES TO 1 AT 20 DEGREES (ABSOLUTE VALUE OF LATITUDE)
C


      TWSTD  = FIVE

      ROOT2I = ONE/SQRT(TWO)



C
C     NOTE: MAKING BC=BF=1/SQRT(2) IN STRIP 10 DEG EITHER SIDE OF EQU:
C
C
C     THIS IS AN EFFORT TO AVOID PURELY DIVERGENT ANALYSIS INCREMENTS
C     AT THE EQUATOR -- SHOULD MAKE ANALYSIS MORE "SCM-LIKE" IN THAT
C     REGION (I.E. EQUALLY ROTATIONAL AND DIVERGENT COMPONENTS TO THE
C     ANALYSIS INCREMENTS)
C

      DO 30 I=1,NWF

              IF(ABS(WLAT(I)) .LE. TEN) THEN

                        WSTD(I) = TWSTD

                        BC(I)   = ROOT2I
                        BF(I)   = ROOT2I

              ELSE IF(ABS(WLAT(I)) .LE. TWENTY) THEN

                        A       = TENTH*ABS( WLAT(I) ) - ONE
                        WSTD(I) = (ONE-A) * TWSTD + A * WSTD(I)

              END IF

 30   CONTINUE


      RETURN
      END
