!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 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