C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE PRODQ2(LMHK,DTQ2,USTAR,GM,GH,EL,Q2)
C ******************************************************************
C * *
C * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
C * *
C ******************************************************************
C-----------------------------------------------------------------------
INCLUDE "parmeta"
# include "sp.h"
C-----------------------------------------------------------------------
P A R A M E T E R
&(LM1=LM-1)
C-----------------------------------------------------------------------
P A R A M E T E R
&(EPSQ2=0.2,EPSL=0.32,EPSTRB=1.E-24,EPS1=1.E-12,EPS2=0.)
C-----------------------------------------------------------------------
P A R A M E T E R
C-----------------------------------------------------------------------
&(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)
C-----------------------------------------------------------------------
CUN &(G=9.8,BETA=1./270.,BTG=BETA*G
CUN &,PRT=1.0,GAM1=0.2222222222222222222
CUN &,A1=0.3310949523016403346,A2=0.8273378704055731278
CUN &,B1=5.959709141429526024,B2=3.626088092074591135
CUN &,C1=-0.3330651924968952113)
C-----------------------------------------------------------------------
CMY &(G=9.8,BETA=1./270.,BTG=BETA*G
CMY &,PRT=0.8,GAM1=0.2222222222222222222
CMY &,A1=0.9222222350809054114,A2=0.7350190142719400952
CMY &,B1=16.60000023145629741,B2=10.10000014082581951
CMY &,C1=0.0805318118080613468)
C-----------------------------------------------------------------------
P A R A M E T E R
&(RB1=1./B1
C--------------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
C--------------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
C--------------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
C--------------FORBIDDEN TURBULENCE AREA--------------------------------
&,REQU=-AEQH/AEQM*1.02,EPSGH=1.E-9)
C-----------------------------------------------------------------------
D I M E N S I O N
& Q2 (LM)
D I M E N S I O N
& GM (LM1),GH (LM1),EL (LM1)
C-----------------------------------------------------------------------
C***********************************************************************
LMHM=LMHK-1
C
DO 150 L=1,LMHM
GML=GM(L)
GHL=GH(L)
C--------------COEFFICIENTS OF THE EQUILIBRIUM EQUATION-----------------
AEQU=(AEQM*GML+AEQH*GHL)*GHL
BEQU= BEQM*GML+BEQH*GHL
C--------------EQUILIBRIUM SOLUTION FOR L/Q-----------------------------
EQOL2=-0.5*BEQU+SQRT(BEQU*BEQU*0.25-AEQU)
C--------------IS THERE PRODUCTION/DISSIPATION ?------------------------
IF((GML+GHL*GHL.LE.EPSTRB )
& .OR.(GHL.GE.EPSGH.AND.GML/GHL.LE.REQU)
& .OR.(EQOL2.LE.EPS2) ) THEN
C--------------NO TURBULENCE--------------------------------------------
Q2(L)=EPSQ2
EL(L)=EPSL
C--------------END OF THE NO TURBULENCE BRANCH--------------------------
ELSE
C--------------COEFFICIENTS OF THE TERMS IN THE NUMERATOR---------------
ANUM=(ANMM*GML+ANMH*GHL)*GHL
BNUM= BNMM*GML+BNMH*GHL
C--------------COEFFICIENTS OF THE TERMS IN THE DENOMINATOR-------------
ADEN=(ADNM*GML+ADNH*GHL)*GHL
BDEN= BDNM*GML+BDNH*GHL
CDEN= 1.
C--------------COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.------
ARHS=-(ANUM*BDEN-BNUM*ADEN)*2.
BRHS=- ANUM*4.
CRHS=- BNUM*2.
C--------------INITIAL VALUE OF L/Q-------------------------------------
DLOQ1=EL(L)/SQRT(Q2(L))
C--------------FIRST ITERATION FOR L/Q, RHS=0---------------------------
ELOQ21=1./EQOL2
ELOQ11=SQRT(ELOQ21)
ELOQ31=ELOQ21*ELOQ11
ELOQ41=ELOQ21*ELOQ21
ELOQ51=ELOQ21*ELOQ31
C--------------1./DENOMINATOR-------------------------------------------
RDEN1=1./(ADEN*ELOQ41+BDEN*ELOQ21+CDEN)
C--------------D(RHS)/D(L/Q)--------------------------------------------
RHSP1= (ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1
C--------------FIRST-GUESS SOLUTION-------------------------------------
ELOQ12=ELOQ11+(DLOQ1-ELOQ11)*EXP(RHSP1*DTQ2)
C-----------------------------------------------------------------------
ELOQ12=AMAX1(ELOQ12,EPS1)
C--------------SECOND ITERATION FOR L/Q---------------------------------
ELOQ22=ELOQ12*ELOQ12
ELOQ32=ELOQ22*ELOQ12
ELOQ42=ELOQ22*ELOQ22
ELOQ52=ELOQ22*ELOQ32
C--------------1./DENOMINATOR-------------------------------------------
RDEN2=1./(ADEN*ELOQ42+BDEN*ELOQ22+CDEN)
C-----------------------------------------------------------------------
RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1
RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2
RHST2=RHS2/RHSP2
C--------------CORRECTED SOLUTION---------------------------------------
ELOQ13=ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*EXP(RHSP2*DTQ2)
C-----------------------------------------------------------------------
ELOQ13=AMAX1(ELOQ13,EPS1)
C--------------TWO ITERATIONS IS ENOUGH IN MOST CASES ...---------------
ELOQN=ELOQ13
C-----------------------------------------------------------------------
IF(ELOQN.GT.EPS1)THEN
Q2(L)=EL(L)*EL(L)/(ELOQN*ELOQN)
Q2(L)=AMAX1(Q2(L),EPSQ2)
ELSE
Q2(L)=EPSQ2
ENDIF
C--------------END OF TURBULENT BRANCH----------------------------------
ENDIF
C--------------END OF PRODUCTION/DISSIPATION LOOP-----------------------
150 CONTINUE
C--------------LOWER BOUNDARY CONDITION FOR Q2--------------------------
Q2(LMHK)=AMAX1(B1**(2./3.)*USTAR*USTAR,EPSQ2)
C-----------------------------------------------------------------------
RETURN
END