#include "wwm_functions.h"
      SUBROUTINE TAUHF_WAM(ML)

! ----------------------------------------------------------------------

!**** *TAUHF* - COMPUTATION OF HIGH-FREQUENCY STRESS.

!     PETER A.E.M. JANSSEN    KNMI      OCTOBER 90

!*    PURPOSE.
!     ---------

!       COMPUTE HIGH-FREQUENCY WAVE STRESS
!       FOR BOTH ECMWF PHYSICS AND METREO FRANCE PHYSICS.

!**   INTERFACE.
!     ----------

!       *CALL* *TAUHF(ML)*
!             *ML*  NUMBER OF FREQUENCIES.

!     METHOD.
!     -------

!       SEE REFERENCE FOR WAVE STRESS CALCULATION.

!     EXTERNALS.
!     ----------

!       NONE.

!     REFERENCE.
!     ----------

!       FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.

! ----------------------------------------------------------------------

!      USE YOWFRED  , ONLY : FR
!      USE YOWCOUP  , ONLY : BETAMAX  ,ZALP     ,ALPHA    ,XKAPPA
!      USE YOWPCONS , ONLY : G        ,ZPI
!      USE YOWTABL  , ONLY : IUSTAR   ,IALPHA   ,USTARM   ,TAUHFT   ,
!     &            DELUST   ,DELALP

       USE DATAPOOL, ONLY : FR, WETAIL, FRTAIL, WP1TAIL, ISHALLO, ILEVTAIL, &
     &                      DFIM, DFIMOFR, DFFR, DFFR2, WK, DELTAIL, IPHYS, &
     &                      IUSTAR, IALPHA, USTARM, TAUHFT, STAT, TAUHFT2, &
     &                      DELUST, DELALP, ALPHA, BETAMAX, RKIND, JTOT, &
     &                      XKAPPA, ZALP, LOUTWAM, &
     &                      DELTH => DDIR, &
     &                      G => G9, &
     &                      ZPI => PI2, &
     &                      EPSMIN => SMALL, &
     &                      NANG => MDC, &
     &                      NFRE => MSC, &
     &                      INDEP => DEP, &
     &                      ZERO, ONE, &
     &                      SRCDBG
#ifdef MPI_PARALL_GRID
      USE DATAPOOL, ONLY : COMM, MYRANK
#endif


      IMPLICIT NONE

! ----------------------------------------------------------------------

      INTEGER, INTENT(IN)  :: ML
      INTEGER              :: I, J, K, L, M

!      REAL(rkind), ALLOCATABLE :: W(:)
      REAL(rkind) :: ALPHAM, ALPHAMCOEF, CONST1, OMEGAC, X0, UST, Z0, OMEGACC, YC
      REAL(rkind) :: DELY, OMEGA, CM, ZX, ZARG, ZMU, ZLOG, Y, ZBETA
      REAL(rkind) :: UST0, XLEVTAIL, TAUW0, TAUW, W(JTOT)
      integer istat

! ----------------------------------------------------------------------

!*    1. PRELIMINARY CALCULATIONS.
!        -------------------------

      ALPHAMCOEF = 40.
      ALPHAM = ALPHAMCOEF*ALPHA
      DELUST = USTARM/REAL(IUSTAR)
      DELALP = ALPHAM/REAL(IALPHA)
      DELTAIL= ALPHAM/REAL(ILEVTAIL)

      CONST1 = BETAMAX/XKAPPA**2

!      GOTO 100

      !ALLOCATE(W(JTOT))
      W=1.
      W(1)=0.5
      W(JTOT)=0.5

      IF (LOUTWAM) WRITE(111111,'(A20)') 'TAUHF'

      IF (LOUTWAM) WRITE(111111,'(10F15.10)') ALPHAM, DELUST, DELALP, DELTAIL, CONST1, SUM(W)


!     TABLES FOR ECMWF PHYSICS:
!     -------------------------

!      IF(.NOT.ALLOCATED(TAUHFT)) ALLOCATE(TAUHFT(0:IUSTAR,0:IALPHA,ML))

      IF (LOUTWAM) WRITE(111111,'(A20,I10)') 'TOTAL NUMBER OF ENTRIES', ML*IALPHA*IUSTAR

!!$OMP PARALLEL DEFAULT(NONE) &
!!$OMP&         SHARED(ML,FR) &
!!$OMP&         PRIVATE(M,L,K,J,OMEGAC,X0,UST,Z0,OMEGACC,&
!!$OMP&         YC,DELY,Y,OMEGA,CM,ZX,ZARG,ZMU,ZLOG,ZBETA,&
!!$OMP&         W,DELALP,ZALP,TAUHFT,DELUST,CONST1,ALPHA)
!!$OMP DO SCHEDULE (DYNAMIC)
      DO M=1,ML

        OMEGAC = ZPI*FR(M)

        DO L=0,IALPHA
          DO K=0,IUSTAR
            TAUHFT(K,L,M) = 0.
          ENDDO
        ENDDO

!*    2. CALCULATE HIGH-FREQUENCY CONTRIBUTION TO STRESS.
!        ------------------------------------------------

        X0 = 0.05
        DO L=0,IALPHA
          DO K=0,IUSTAR
            UST      = MAX(REAL(K)*DELUST,0.000001)
            Z0       = UST**2*(ALPHA+REAL(L)*DELALP)/G
            OMEGACC  = MAX(OMEGAC,X0*G/UST)
            YC       = OMEGACC*SQRT(Z0/G)
            DELY     = MAX((1.-YC)/REAL(JTOT-1),0.)
            DO J=1,JTOT
              Y        = YC+REAL(J-1)*DELY
              OMEGA    = Y*SQRT(G/Z0)
              CM       = G/OMEGA
              ZX       = UST/CM +ZALP
              ZARG     = MIN(XKAPPA/ZX,20.)
              ZMU      = MIN(G*Z0/CM**2*EXP(ZARG),1.)

              ZLOG         = MIN(LOG(ZMU),0.)
              ZBETA        = CONST1*ZMU*ZLOG**4
              TAUHFT(K,L,M)= TAUHFT(K,L,M)+W(J)*ZBETA/Y*DELY
            ENDDO
!            WRITE(111111,'(3I10,10F15.7)') L,K,J,DELY,ZMU,TAUHFT(K,L,M)
          ENDDO
        ENDDO

      ENDDO
!!$OMP END PARALLEL



!     TABLES FOR METEO FRANCE PHYSICS:
!     -------------------------------

!      IF(.NOT.ALLOCATED(TAUHFT2)) &
!     & ALLOCATE(TAUHFT2(0:IUSTAR,0:IALPHA,0:ILEVTAIL))


!*    3. CALCULATE HIGH-FREQUENCY CONTRIBUTION TO STRESS.
!        ------------------------------------------------

!!$OMP PARALLEL DEFAULT(NONE) &
!!$OMP&         SHARED(FR,ML) &
!!$OMP&         PRIVATE(L,K,I,J,OMEGAC,X0,UST0,UST,Z0,OMEGACC,&
!!$OMP&         YC,DELY,XLEVTAIL,TAUW0,TAUW,Y,OMEGA,CM,ZX,ZARG,&
!!$OMP&         ZMU,ZLOG,ZBETA,CONST1,W,ZALP,DELALP,DELTAIL,&
!!$OMP&         DELUST,ALPHA,TAUHFT2)
!!$OMP DO SCHEDULE (DYNAMIC)
      DO L=0,IALPHA
        OMEGAC = ZPI*FR(ML)
        X0 = 0.05
        DO K=0,IUSTAR
          UST0     = MAX(REAL(K)*DELUST,0.000001)
          UST      = UST0
          Z0       = UST0**2*(ALPHA+FLOAT(L)*DELALP)/G
          OMEGACC  = MAX(OMEGAC,X0*G/UST)
          YC       = OMEGACC*SQRT(Z0/G)
          DELY     = MAX((1.-YC)/REAL(JTOT),0.)
          DO I=0,ILEVTAIL
            TAUHFT2(K,L,I)=0.
            XLEVTAIL=REAL(I)*DELTAIL
            TAUW0=UST0**2
            TAUW=TAUW0
            DO J=1,JTOT
              Y        = YC+REAL(J-1)*DELY
              OMEGA    = Y*SQRT(G/Z0)
              CM       = G/OMEGA
              ZX       = UST/CM +ZALP
              ZARG     = MIN(XKAPPA/ZX,20.)
              ZMU      = MIN(G*Z0/CM**2*EXP(ZARG),1.)

              ZLOG         = MIN(LOG(ZMU),0.)
              ZBETA        = CONST1*ZMU*ZLOG**4
              TAUHFT2(K,L,I)= &
     &           TAUHFT2(K,L,I)+W(J)*ZBETA*(UST/UST0)**2/Y*DELY 
              TAUW         =TAUW-W(J)*UST**2*ZBETA*XLEVTAIL/Y*DELY
              UST          =SQRT(MAX(TAUW,0.))
            ENDDO
          ENDDO
        ENDDO
      ENDDO
!!$OMP END PARALLEL

      IF (LOUTWAM) WRITE(111111,'(A10,I10)') 'IPHYS=', IPHYS

      IF (IPHYS == 0) THEN
#ifdef MPI_PARALL_GRID
        if (myrank == 0 ) then
          WRITE(5011) DELALP, DELUST, DELTAIL
          WRITE(5011) TAUHFT
        endif
#endif
        IF (LOUTWAM) WRITE(111111,'(F20.10)') SUM(TAUHFT)
      ELSE
#ifdef MPI_PARALL_GRID
        if (myrank ==0 ) then
          WRITE(5011) DELALP, DELUST, DELTAIL
          WRITE(5011) TAUHFT, TAUHFT2, TAUW
        endif
#endif
        IF (LOUTWAM) WRITE(111111,'(3F20.10)') DELTAIL, SUM(TAUHFT), SUM(TAUHFT2) 
      ENDIF

!      DO M=1,NFRE
!        WRITE(100000,*) M, FR(M)
!      ENDDO

!      DO M=1,NFRE
!        DO L=0,IALPHA
!          DO K=0,IUSTAR
!            WRITE(100002,*) K,L,M,TAUHFT(K,L,M)
!          ENDDO
!        ENDDO
!      ENDDO

!      DEALLOCATE(W)

      END SUBROUTINE TAUHF_WAM