!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE PRODQ2(LMHK,DTQ2,USTAR,GM,GH,EL,Q2)
!     ******************************************************************
!     *                                                                *
!     *  LEVEL 2.5 Q2 PRODUCTION/DISSIPATION                           *
!     *                                                                *
!     ******************************************************************
!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
#include "sp.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (LM1=LM-1)
!-----------------------------------------------------------------------
    PARAMETER &
    (EPSQ2=0.2,EPSL=0.32,EPSTRB=1.E-24,EPS1=1.E-12,EPS2=0.)
!-----------------------------------------------------------------------
    PARAMETER &
!-----------------------------------------------------------------------
    (G=9.8,BETA=1./270.,BTG=BETA*G &
    ,PRT=1.0,GAM1=0.2222222222222222222 &
    ,A1=0.659888514560862645,A2=0.6574209922667784586 &
    ,B1=11.87799326209552761,B2=7.226971804046074028 &
    ,C1=0.000830955950095854396)
!-----------------------------------------------------------------------
! N     &(G=9.8,BETA=1./270.,BTG=BETA*G
! N     &,PRT=1.0,GAM1=0.2222222222222222222
! N     &,A1=0.3310949523016403346,A2=0.8273378704055731278
! N     &,B1=5.959709141429526024,B2=3.626088092074591135
! N     &,C1=-0.3330651924968952113)
!-----------------------------------------------------------------------
! Y     &(G=9.8,BETA=1./270.,BTG=BETA*G
! Y     &,PRT=0.8,GAM1=0.2222222222222222222
! Y     &,A1=0.9222222350809054114,A2=0.7350190142719400952
! Y     &,B1=16.60000023145629741,B2=10.10000014082581951
! Y     &,C1=0.0805318118080613468)
!-----------------------------------------------------------------------
    PARAMETER &
    (RB1=1./B1 &
!--------------COEFFICIENTS OF THE TERMS IN THE NUMERATOR---------------
    ,ANMM=-3.*A1*A2*(3.*A2+3.*B2*C1+18.*A1*C1-B2)*BTG &
    ,ANMH=-9.*A1*A2*A2*BTG*BTG &
    ,BNMM=    A1*(1.-3.*C1) &
    ,BNMH=   -A2*BTG &
!--------------COEFFICIENTS OF THE TERMS IN THE DENOMINATOR-------------
    ,ADNM=18.*A1*A1*A2*(B2-3.*A2)*BTG &
    ,ADNH= 9.*A1*A2*A2*(12.*A1+3.*B2)*BTG*BTG &
    ,BDNM= 6.*A1*A1 &
    ,BDNH= 3.*A2*(7.*A1+B2)*BTG &
!--------------COEFFICIENTS OF THE EQUILIBRIUM EQUATION-----------------
    ,AEQM= 3.*A1*A2*B1*(3.*A2+3.*B2*C1+18.*A1*C1-B2)*BTG &
    +18.*A1*A1*A2*(B2-3.*A2)*BTG &
    ,AEQH= 9.*A1*A2*A2*B1*BTG*BTG+9.*A1*A2*A2*(12.*A1+3.*B2)*BTG*BTG &
    ,BEQM=-A1*B1*(1.-3.*C1)+6.*A1*A1 &
    ,BEQH= A2*B1*BTG+3.*A2*(7.*A1+B2)*BTG &
!--------------FORBIDDEN TURBULENCE AREA--------------------------------
    ,REQU=-AEQH/AEQM*1.02,EPSGH=1.E-9)
!-----------------------------------------------------------------------
    DIMENSION &
    Q2    (LM)
    DIMENSION &
    GM    (LM1),GH    (LM1),EL    (LM1)
!-----------------------------------------------------------------------
!***********************************************************************
    LMHM=LMHK-1

    DO 150 L=1,LMHM
        GML=GM(L)
        GHL=GH(L)
    !--------------COEFFICIENTS OF THE EQUILIBRIUM EQUATION-----------------
        AEQU=(AEQM*GML+AEQH*GHL)*GHL
        BEQU= BEQM*GML+BEQH*GHL
    !--------------EQUILIBRIUM SOLUTION FOR L/Q-----------------------------
        EQOL2=-0.5*BEQU+SQRT(BEQU*BEQU*0.25-AEQU)
    !--------------IS THERE PRODUCTION/DISSIPATION ?------------------------
        IF((GML+GHL*GHL <= EPSTRB           ) &
         .OR. (GHL >= EPSGH .AND. GML/GHL <= REQU) &
         .OR. (EQOL2 <= EPS2)                  )    THEN
        !--------------NO TURBULENCE--------------------------------------------
            Q2(L)=EPSQ2
            EL(L)=EPSL
        !--------------END OF THE NO TURBULENCE BRANCH--------------------------
        ELSE
        !--------------COEFFICIENTS OF THE TERMS IN THE NUMERATOR---------------
            ANUM=(ANMM*GML+ANMH*GHL)*GHL
            BNUM= BNMM*GML+BNMH*GHL
        !--------------COEFFICIENTS OF THE TERMS IN THE DENOMINATOR-------------
            ADEN=(ADNM*GML+ADNH*GHL)*GHL
            BDEN= BDNM*GML+BDNH*GHL
            CDEN= 1.
        !--------------COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.------
            ARHS=-(ANUM*BDEN-BNUM*ADEN)*2.
            BRHS=- ANUM*4.
            CRHS=- BNUM*2.
        !--------------INITIAL VALUE OF L/Q-------------------------------------
            DLOQ1=EL(L)/SQRT(Q2(L))
        !--------------FIRST ITERATION FOR L/Q, RHS=0---------------------------
            ELOQ21=1./EQOL2
            ELOQ11=SQRT(ELOQ21)
            ELOQ31=ELOQ21*ELOQ11
            ELOQ41=ELOQ21*ELOQ21
            ELOQ51=ELOQ21*ELOQ31
        !--------------1./DENOMINATOR-------------------------------------------
            RDEN1=1./(ADEN*ELOQ41+BDEN*ELOQ21+CDEN)
        !--------------D(RHS)/D(L/Q)--------------------------------------------
            RHSP1= (ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1
        !--------------FIRST-GUESS SOLUTION-------------------------------------
            ELOQ12=ELOQ11+(DLOQ1-ELOQ11)*EXP(RHSP1*DTQ2)
        !-----------------------------------------------------------------------
            ELOQ12=AMAX1(ELOQ12,EPS1)
        !--------------SECOND ITERATION FOR L/Q---------------------------------
            ELOQ22=ELOQ12*ELOQ12
            ELOQ32=ELOQ22*ELOQ12
            ELOQ42=ELOQ22*ELOQ22
            ELOQ52=ELOQ22*ELOQ32
        !--------------1./DENOMINATOR-------------------------------------------
            RDEN2=1./(ADEN*ELOQ42+BDEN*ELOQ22+CDEN)
        !-----------------------------------------------------------------------
            RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1
            RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2
            RHST2=RHS2/RHSP2
        !--------------CORRECTED SOLUTION---------------------------------------
            ELOQ13=ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*EXP(RHSP2*DTQ2)
        !-----------------------------------------------------------------------
            ELOQ13=AMAX1(ELOQ13,EPS1)
        !--------------TWO ITERATIONS IS ENOUGH IN MOST CASES ...---------------
            ELOQN=ELOQ13
        !-----------------------------------------------------------------------
            IF(ELOQN > EPS1)THEN
                Q2(L)=EL(L)*EL(L)/(ELOQN*ELOQN)
                Q2(L)=AMAX1(Q2(L),EPSQ2)
            ELSE
                Q2(L)=EPSQ2
            ENDIF
        !--------------END OF TURBULENT BRANCH----------------------------------
        ENDIF
    !--------------END OF PRODUCTION/DISSIPATION LOOP-----------------------
    150 END DO
!--------------LOWER BOUNDARY CONDITION FOR Q2--------------------------
    Q2(LMHK)=AMAX1(B1**(2./3.)*USTAR*USTAR,EPSQ2)
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE PRODQ2
