CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of WVARZ from wvarz.f
CMAX: ----------------------------------------------------
      SUBROUTINE wvarz(ireg,nwf,siny,cosy,wlat,wstd,wprs,bc,bf,scale)
      include '/usr/include/cm/CMF_defs.h'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C         Optimum Interpolation Analysis Version 1.1
C
C
C WVARZ FORMS THE ARRAY OF WIND FCST ERROR VARIANCES (IN CORRECT UNITS
C (UPPER-AIR ANALYSIS)
C
C  INPUT ARGUMENTS:
C
C     NWF      = NW + NH*NV
C     WSTD     = HEIGHT VARIANCES AT GRID POINTS
C
C  OUTPUT ARGUMENTS:
C
C                                  ABS(LAT)
C     BC       = SGN * ( 1 - EXP (--------- ) , ABS(LAT) < 25
C                                    25
C
C                WHERE SGN=+1 (NH) AND SGN =-1 (SH)
C
C
C                   SINY(I)
C     BF       =  -----------     ,  HEMISPHERE SIGN FACTOR
C                 ABS(SINY(I)
C
C     SCALE    = DIMENSIONAL FACTOR FROM CORRELATION TABLE
C     WSTD     = WIND ERROR VARIANCE (IN UNITS OF M SEC-1)
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     FOR CURRENT TESTS...USING TWSTD = 5 M SEC-1
C
C
C     NEW ROUTINE FOR O9X OI, DEVELOPED BY S.BLOOM   4/90
C     Full implementation to version O90 (1.0)       9/90
C
C     Modifications for OI 1.1:
C
C     Changes to allow for > 12 levels. (S. Bloom)   3/91
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PARAMETER (maxplvl = 17, maxtwerr = 18)
      DIMENSION wstd(:),wlat(:),siny(:),cosy(:),bf(:),bc(:)
      DIMENSION twstd(700),plvl(maxplvl),twerr(maxtwerr),wprs(:)
      CHARACTER * 60 message
      COMMON / GRIDXX / IM, JNP, MLEV
      COMMON / PRINDEX / PINDX (1050)
      COMMON / CORTABS / HCORFE (1000, 3, 20, 2)
!hpf$ distribute bc(block)
!hpf$ distribute bf(block)
!hpf$ distribute cosy(block)
!hpf$ distribute hcorfe(block,block,block,block)
!hpf$ distribute pindx(block)
!hpf$ distribute plvl(block)
!hpf$ distribute siny(block)
!hpf$ distribute twerr(block)
!hpf$ distribute twstd(block)
!hpf$ distribute wlat(block)
!hpf$ distribute wprs(block)
!hpf$ distribute wstd(block)
      REAL  DEN100(nwf)
!hpf$ distribute DEN100(block)
      REAL  A100(nwf)
!hpf$ distribute A100(block)
      REAL  BCI100(nwf)
!hpf$ distribute BCI100(block)
      DATA ZERO / 0.0 /
      DATA ONE / 1.0 /
      DATA ONEPT5 / 1.5 /
      DATA TWO / 2.0 /
      DATA TEN / 10.0 /
      DATA ONE5 / 15.0 /
      DATA TWENTY / 20.0 /
      DATA TWNTY5 / 25.0 /
      DATA TENTH / 0.1 /
      DATA TENM3 / 0.001 /
      DATA TENM5 / 0.00001 /
c     DATA PLVL/ 60.e0, 85.e0,125.e0,175.e0,225.e0,275.e0,350.e0,
c     DATA PLVL/ 60.e0, 85.e0,125.e0,175.e0,225.e0,275.e0,350.e0,
      DATA PLVL / 0.7e0, 1.5e0, 3.5e0, 7.5e0, 20.e0, 40.e0, 60.e0, 85.e0
     &, 125.e0, 175.e0, 225.e0, 275.e0, 350.e0, 450.e0, 600.e0, 775.e0, 
     &925.e0 /
      DATA TWERR / 8.0e0, 7.3e0, 7.3e0, 7.2e0, 7.1e0, 7.1e0, 8.0e0, 7.3e
     &0, 7.3e0, 7.3e0, 7.4e0, 7.0e0, 6.2e0, 5.2e0, 4.5e0, 4.4e0, 4.4e0, 
     &4.3e0 /
C**************************************************************
C****                                                    ******
C****   THIS IS A CRUDE FIRST SHOT AT ESTIMATING TWERR   ******
C****   FOR STRATOSPHERIC ANALYSIS -- IT SHOULD BE       ******
C****   CHECKED FOR CONSISTENCY DURING ANALYSIS TESTS    ******
C****                                                    ******
C****   (NB I'VE ASSUMED MAXIMA NEAR TROPOPAUSE AND TOP) ******
C****                                                    ******
C**************************************************************
c     DATA PLVL/  0.7e0, 1.5e0, 3.5e0, 7.5e0, 20.e0, 40.e0,
c    .            60.e0, 85.e0,125.e0,175.e0,225.e0,275.e0,
c    .           350.e0,450.e0,600.e0,775.e0,925.e0/
c     DATA TWERR/8.0e0,7.3e0,7.3e0,7.2e0,7.1e0,7.1e0,
c    .           7.1e0,7.2e0,7.3e0,7.3e0,7.4e0,7.0e0,
c    .           6.2e0,5.2e0,4.5e0,4.4e0,4.4e0,4.3e0/
      CALL getconh('GRADEI  ',gradei)
      CALL getconh('PI180   ',pi180)
      CALL getconh('OMEGA2  ',omega2)
      CALL getconh('RADEI   ',radei)
      icor = 2
      message(1:60) = ' '
C
C.... OBTAIN SCALE FACTOR FROM CORRELATION TABLE
C
      scale = hcorfe(1,3,ireg,icor)
      IF (scale .GE. zero) THEN
        message = ' WVARZ: Bad scale value detected. '
        CALL flash(5,message,1)
      ENDIF
      scale = sqrt(-scale)
C
C.... HERE WE CONVERT WSTD INTO A GEOSTROPHIC ESTIMATE OF WIND
C     ERROR VARIANCE...USING FORMULA:
C
C
C
C        WSTD(GEOST) = SCALE * G * WSTD(OLD) / (RADE * ABS(F) )
C
C     NOTE, THE UNITS OF WSTD(NEW) ARE M SEC-1
C     AND THAT WE INITIALLY PUT GEOSTROPHIC BALANCE ESTIMATES IN
C     WSTD FOR ALL LATITUDES
C
C     ALSO, GEOSTROPHIC ESTIMATES USE CORIOLIS PARAMETER 'PEGGED'
C     TO ABSOLUTE LATITUDES NO SMALLER THAN 25 DEGREES.
C
      coef = gradei * scale / omega2
      cfsin = one / sin(pi180 * twnty5)
      den100(1:nwf) = one / abs(siny(:nwf))
      bf(:nwf) = siny(:nwf) * den100(1:nwf)
      WHERE (abs(wlat(:nwf)) .LE. twnty5)
        den100(1:nwf) = cfsin
      ENDWHERE
      wstd(:nwf) = coef * wstd(:nwf) * den100(1:nwf)
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
C*************************************************************
C*************************************************************
C
C     TWO CHOICES TO OBTAIN TWSTD
C
C             1.  ORIGINAL WAY IN OLD ZUVCOV - VERY 'KLUDGY'
C             2.  'IMPROVED' METHOD USING APPROPRIATE VERTICAL
C                 INTERPOLATION
C
C*************************************************************
C*************************************************************
C
C     --------
C     METHOD 1
C     --------
C
      FORALL (i = 1:nwf, wprs(i) .GE. plvl(maxplvl)) twstd(i) 
     &    = twerr(maxtwerr)
C         IF( WPRS(I) .LT. PLVL(1) )  TWSTD(I) = TWERR(1)
      FORALL (i = 1:nwf, wprs(i) .LT. plvl(maxplvl - mlev + 2)) twstd(i)
     &    =twerr(maxplvl-mlev+2)
      DO i = 1 , nwf
        DO jkl = 2 , maxplvl
          kkm = jkl - 1
          IF (wprs(i) .GE. plvl(kkm) .AND. wprs(i) .LT. plvl(jkl)) THEN
            twstd(i) = twerr(jkl)
          ENDIF
        ENDDO
      ENDDO
C
C     --------
C     METHOD 2    (USING OBERTRP METHOD)
C     --------
C
C     DO 28 I=1,NWF
C
C         hl    = pindx2( wprs(i) )
C         LEVUP = hl + TENM3
C         LEVUP = MAX0(LEVUP,1)
C         LEVDN = MIN0(LEVUP+1,MLEV)
C
C         TWSTD(I) = TWERR(LEVUP) + (hl- LEVUP) *
C    1               ( TWERR(LEVDN) - TWERR(LEVUP) )
C
C28   CONTINUE
      WHERE (abs(wlat(:nwf)) .LE. one5)
        wstd(:nwf) = twstd(:nwf)
      ENDWHERE
      WHERE (.NOT.abs(wlat(:nwf)) .LE. one5 .AND. abs(wlat(:nwf)) 
     &      .LE. twnty5)
        a100(1:nwf) = tenth * abs(wlat(:nwf)) - onept5
        wstd(:nwf) = (one - a100(1:nwf)) * twstd(:nwf) + a100(1:nwf) 
     &      * wstd(:nwf)
      ENDWHERE
C
C     LOAD WIND-HEIGHT DECOUPLING FACTOR INTO BC (NEEDED IN WMCOV AND
C                                                           WWCOV)
C
      scexp = -one / twnty5
      bci100(1:nwf) = one - exp(scexp * abs(wlat(:nwf)))
      WHERE (bci100(1:nwf) .EQ. zero)
        bci100(1:nwf) = tenm5
      ENDWHERE
      bc(:nwf) = bci100(1:nwf) * bf(:nwf)
      WHERE (abs(wlat(:nwf)) .LT. 15.)
        bf(:nwf) = 0.
      ENDWHERE
      RETURN
      END
