SUBROUTINE CALTAU(TAUX,TAUY)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    CALTAU      COMPUTE U AND V WIND STRESSES
C   PRGRMMR: TREADON         ORG: W/NP2      DATE: 93-09-01
C     
C ABSTRACT:  THIS ROUTINE COMPUTES SURFACE LAYER U AND V
C   WIND COMPONENT STRESSES USING K THEORY AS PRESENTED
C   IN SECTION 8.4 OF "NUMBERICAL PREDICTION AND DYNAMIC
C   METEOROLOGY" BY HALTINER AND WILLIAMS (1980, JOHN WILEY
C   & SONS).
C   .     
C     
C PROGRAM HISTORY LOG:
C   93-09-01  RUSS TREADON
C   98-06-11  T BLACK - CONVERSION FROM 1-D TO 2-D
C   00-01-04  JIM TUCCILLO - MPI VERSION
C     
C USAGE:    CALL CALTAU(TAUX,TAUY)
C   INPUT ARGUMENT LIST:
C     NONE     
C
C   OUTPUT ARGUMENT LIST: 
C     TAUX     - SUFACE LAYER U COMPONENT WIND STRESS.
C     TAUY     - SUFACE LAYER V COMPONENT WIND STRESS.
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       CLMAX
C       MIXLEN
C       H2V
C
C     LIBRARY:
C       COMMON   - EXTRA
C                  LOOPS
C                  VRBLS
C                  PVRBLS
C                  MAPOT
C                  OPTIONS
C                  INDX
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C
C     INCLUDE/SET PARAMETERS.
C     
      INCLUDE "parmeta"
      INCLUDE "parmout"
      INCLUDE "params"
C     
C     DECLARE VARIABLES.
      INTEGER KK(4)
      REAL TAUX(IM,JM),TAUY(IM,JM)
      REAL EL0(IM,JM),EL(IM,JM,LM)
      REAL EGRIDU(IM,JM),EGRIDV(IM,JM),EGRID4(IM,JM),EGRID5(IM,JM)
C     
C     INCLUDE COMMON BLOCKS.
C     
      INCLUDE "EXTRA.comm"
      INCLUDE "LOOPS.comm"
      INCLUDE "VRBLS.comm"
      INCLUDE "PVRBLS.comm"
      INCLUDE "MAPOT.comm"
      INCLUDE "MASKS.comm"
      INCLUDE "OPTIONS.comm"
      INCLUDE "INDX.comm"
      INCLUDE "CTLBLK.comm"
C     
C     
C********************************************************************
C     START CALTAU HERE.
C     
C     COMPUTE MASTER LENGTH SCALE.
C
      CALL CLMAX(DETA,PDSL,HTM,Q2,ZINT,SM,ZINT(1,1,LP1),
     1     LMH,IM,JM,LM,LP1,
     2     EL0,EGRIDU,EGRIDV,EGRID4,EGRID5)
      CALL MIXLEN(ZINT,T,PDSL,AETA,PT,Q2,ZINT(1,1,LP1),HTM
     1,    EL0,LM,LM1,LP1,IM,JM,EL)
C     
C     INITIALIZE OUTPUT AND WORK ARRAY TO ZERO.
C     
      DO J=JSTA,JEND
      DO I=1,IM
        EGRIDU(I,J) = D00
        EGRIDV(I,J) = D00
        TAUX(I,J)   = SPVAL
        TAUY(I,J)   = SPVAL
      ENDDO
      ENDDO
C
      CALL EXCH(UZ0)
      CALL EXCH(VZ0)
C     
C     COMPUTE SURFACE LAYER U AND V WIND STRESSES.
C
C     ASSUME THAT U AND V HAVE UPDATED HALOS
C
      doout30: DO J=JSTA_M,JEND_M
      doin30: DO I=2,IM-1
C
        LMHK = LMH(I,J)
C
C       COMPUTE THICKNESS OF LAYER AT MASS POINT.
C
        DZ  = D50*(ZINT(I,J,LMHK)-ZINT(I,J,LMHK+1))
        DZ  = DZ-Z0(I,J)
        RDZ = 1./DZ
C
C        COMPUTE REPRESENTATIVE AIR DENSITY.
C
        PSFC = PD(I,J)+PT
        TV   = (H1+D608*Q(I,J,LMHK))*T(I,J,LMHK)
        RHO  = PSFC/(RD*TV)
C     
C        COMPUTE A MEAN MASS POINT WIND IN THE 
C        FIRST ATMOSPHERIC ETA LAYER.
C
        LMHK  = LMH(I,J)
        IE=I+IHE(J)
        IW=I+IHW(J)
        SUMU=U(IE,J,LMV(IE,J))+U(IW,J,LMV(IW,J))+U(I,J-1,LMV(I,J-1)) 
     1        +U(I,J+1,LMV(I,J+1))
        SUMV=V(IE,J,LMV(IE,J))+V(IW,J,LMV(IW,J))+V(I,J-1,LMV(I,J-1)) 
     1        +V(I,J+1,LMV(I,J+1))
        ULMH = D25*SUMU
        VLMH = D25*SUMV
C     
C       COMPUTE A MEAN MASS POINT WIND AT HEIGHT Z0.
C
        UZ0H = D25*(UZ0(I,J-1)+UZ0(IW,J)+UZ0(IE,J)+UZ0(I,J+1))
        VZ0H = D25*(VZ0(I,J-1)+VZ0(IW,J)+VZ0(IE,J)+VZ0(I,J+1))
C
C       COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
C
        DELUDZ = (ULMH-UZ0H)*RDZ
        DELVDZ = (VLMH-VZ0H)*RDZ
C     
C       COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
C
        ELSQR       = EL(I,J,LMHK)*EL(I,J,LMHK)
        EGRIDU(I,J) = RHO*ELSQR*DELUDZ*DELUDZ
        EGRIDV(I,J) = RHO*ELSQR*DELVDZ*DELVDZ

CX         WRITE(80,1000) K,LMHK,DZ,RDZ,PSFC,TV
CX         WRITE(80,1000) K,LMHK,ULMH,VLMH,UZ0H,VZ0H
CX         WRITE(80,1000) K,LMHK,RHO,DELUDZ,DELVDZ,EL(I,J,LMHK)
CX         WRITE(80,*)' '
CX 1000    FORMAT(I5,1X,I2,4(G12.6,1X))

C
      END DO doin30
      END DO doout30
C     
C     INTERPOLATE U-V STRESS FROM MASS TO VELOCITY POINTS.
C     
      CALL H2V(EGRIDU,EGRIDV,TAUX,TAUY)
C     
C     END OF ROUTINE.
C     
      RETURN
      END