SUBROUTINE prodq2(lmh,dtq2,ustar,gm,gh,el,q2)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  !  level 2.5 q2 production/dissipation
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: lm1=lm-1
  !
  INTEGER, DIMENSION(im,jm    ), INTENT(IN   ) :: lmh 
  REAL   ,                       INTENT(IN   ) :: dtq2
  REAL   , DIMENSION(im,jm    ), INTENT(IN   ) :: ustar
  REAL   , DIMENSION(im,jm,lm1), INTENT(IN   ) :: gm
  REAL   , DIMENSION(im,jm,lm1), INTENT(IN   ) :: gh
  REAL   , DIMENSION(im,jm,lm1), INTENT(INOUT) :: el
  REAL   , DIMENSION(im,jm,lm ), INTENT(INOUT) :: q2
  !
  REAL, PARAMETER :: epsq2=0.2
  REAL, PARAMETER :: epsl=0.32
  REAL, PARAMETER :: epstrb=1.e-24
  REAL, PARAMETER :: eps1=1.e-12
  REAL, PARAMETER :: eps2=0.
  REAL, PARAMETER :: g=9.8
  REAL, PARAMETER :: beta=1./270.
  REAL, PARAMETER :: btg=beta*g
  REAL, PARAMETER :: prt=1.0
  REAL, PARAMETER :: gam1=0.2222222222222222222
  REAL, PARAMETER :: a1=0.659888514560862645
  REAL, PARAMETER :: a2=0.6574209922667784586
  REAL, PARAMETER :: b1=11.87799326209552761
  REAL, PARAMETER :: b2=7.226971804046074028
  REAL, PARAMETER :: c1=0.000830955950095854396

  REAL, PARAMETER :: rb1=1./b1

  !*** coefficients of the terms in the numerator

  REAL, PARAMETER :: anmm=-3.*a1*a2*(3.*a2+3.*b2*c1+18.*a1*c1-b2)*btg
  REAL, PARAMETER :: anmh=-9.*a1*a2*a2*btg*btg
  REAL, PARAMETER :: bnmm=    a1*(1.-3.*c1)
  REAL, PARAMETER :: bnmh=   -a2*btg

  !*** coefficients of the terms in the denominator

  REAL, PARAMETER :: adnm=18.*a1*a1*a2*(b2-3.*a2)*btg
  REAL, PARAMETER :: adnh= 9.*a1*a2*a2*(12.*a1+3.*b2)*btg*btg
  REAL, PARAMETER :: bdnm= 6.*a1*a1
  REAL, PARAMETER :: bdnh= 3.*a2*(7.*a1+b2)*btg

  !*** coefficients of the equilibrium equation

  REAL, PARAMETER :: aeqm= 3.*a1*a2*b1*(3.*a2+3.*b2*c1+18.*a1*c1-b2)*btg&
       +18.*a1*a1*a2*(b2-3.*a2)*btg
  REAL, PARAMETER :: aeqh= 9.*a1*a2*a2*b1*btg*btg+&
       9.*a1*a2*a2*(12.*a1+3.*b2)*btg*btg
  REAL, PARAMETER :: beqm=-a1*b1*(1.-3.*c1)+6.*a1*a1
  REAL, PARAMETER :: beqh= a2*b1*btg+3.*a2*(7.*a1+b2)*btg

  !*** forbidden turbulence area

  REAL, PARAMETER :: requ=-aeqh/aeqm*1.02
  REAL, PARAMETER :: epsgh=1.e-9

  INTEGER :: i, j, l
  REAL    :: aequ, bequ, eqol2
  REAL    :: anum, bnum, aden, bden, cden
  REAL    :: arhs, brhs, crhs, dloq1
  REAL    :: eloq21, eloq11, eloq31, eloq41, eloq51
  REAL    :: eloq12, eloq22, eloq32, eloq42, eloq52
  REAL    :: eloq13, eloqn
  REAL    :: rden1, rhsp1, rden2, rhs2, rhsp2, rhst2

  DO l=1,lm-1
     DO j = 3, jm-2
        DO i = 1, im
           IF (l <= lmh(i,j)-1) THEN

              !*** coefficients of the equilibrium equation

              aequ=(aeqm*gm(i,j,l)+aeqh*gh(i,j,l))*gh(i,j,l)
              bequ= beqm*gm(i,j,l)+beqh*gh(i,j,l)

              !*** equilibrium solution for l/q

              eqol2=-0.5*bequ+SQRT(bequ*bequ*0.25-aequ)

              !*** is there production/dissipation ?

              IF( (gm(i,j,l)+gh(i,j,l)*gh(i,j,l).LE.epstrb           ) .OR. &
                   (gh(i,j,l).GE.epsgh.AND.gm(i,j,l)/gh(i,j,l).LE.requ) .OR. &
                   (eqol2.LE.eps2)                                    ) THEN

                 !*** no turbulence

                 q2(i,j,l)=epsq2
                 el(i,j,l)=epsl
              ELSE

                 !*** coefficients of the terms in the numerator

                 anum=(anmm*gm(i,j,l)+anmh*gh(i,j,l))*gh(i,j,l)
                 bnum= bnmm*gm(i,j,l)+bnmh*gh(i,j,l)

                 !*** coefficients of the terms in the denominator

                 aden=(adnm*gm(i,j,l)+adnh*gh(i,j,l))*gh(i,j,l)
                 bden= bdnm*gm(i,j,l)+bdnh*gh(i,j,l)
                 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(i,j,l)/SQRT(q2(i,j,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.GT.eps1)THEN
                    q2(i,j,l)=el(i,j,l)*el(i,j,l)/(eloqn*eloqn)
                    q2(i,j,l)=amax1(q2(i,j,l),epsq2)
                 ELSE
                    q2(i,j,l)=epsq2
                 ENDIF
              ENDIF
           END IF
        END DO
     END DO
  END DO

  !*** lower boundary condition for q2

  DO j = 3, jm-2
     DO i = 1, im
        q2(i,j,lmh(i,j))=amax1(b1**(2./3.)*ustar(i,j)*ustar(i,j),epsq2)
     END DO
  END DO
END SUBROUTINE prodq2
