SUBROUTINE PRODQ2(LMHK, DTQ2, USTAR, GM, GH, EL, Q2) !>-------------------------------------------------------------------------------------------------- !> SUBROUTINE PRODQ2 !> !> SUBROUTINE: PRODQ2 - Q2 PRODUCTION/DISSIPATION !> PROGRAMMER: ????? !> ORG: ????? !> DATE: ??-??-?? !> !> ABSTRACT: !> LEVEL 2.5 Q2 PRODUCTION/DISSIPATION !> !> PROGRAM HISTORY LOG: !> 94-??-?? ????? - ORIGINATOR !> 18-01-15 LUCCI - MODERNIZATION OF THE CODE, INCLUDING: !> * F77 TO F90/F95 !> * INDENTATION & UNIFORMIZATION CODE !> * REPLACEMENT OF COMMONS BLOCK FOR MODULES !> * DOCUMENTATION WITH DOXYGEN !> * OPENMP FUNCTIONALITY !> !> INPUT ARGUMENT LIST: !> LMHK - MASS POINT MODEL SURFACE (LOWEST MODEL HEIGHT) !> DTQ2 - PHYSICS TINE STEP !> USTAR - SURFACE FRICTION VELOCITY !> GM - WIND SHEAR !> GH - VERTICAL GRADIENT OS POTENTIAL TEMPERATURE !> !> OUTPUT ARGUMENT LIST: !> EL - !> !> INPUT/OUTPUT ARGUMENT LIST: !> Q2 - TURBULENT KINETIC ENERGY !> !> OUTPUT FILES: !> NONE !> !> USE MODULES: F77KINDS !> PARMETA !> !> DRIVER : TURBL !> !> CALLS : ----- !>-------------------------------------------------------------------------------------------------- USE F77KINDS USE PARMETA ! IMPLICIT NONE ! #include "sp.h" ! INTEGER(KIND=I4KIND), PARAMETER :: LM1 = LM - 1 ! REAL (KIND=R4KIND), PARAMETER :: EPSQ2 = 0.12 REAL (KIND=R4KIND), PARAMETER :: EPSL = 0.32 REAL (KIND=R4KIND), PARAMETER :: EPSTRB = 1.E-24 REAL (KIND=R4KIND), PARAMETER :: EPS1 = 1.E-12 REAL (KIND=R4KIND), PARAMETER :: EPS2 = 0. ! REAL (KIND=R4KIND), PARAMETER :: G = 9.8 REAL (KIND=R4KIND), PARAMETER :: BETA = 1. / 270. REAL (KIND=R4KIND), PARAMETER :: BTG = BETA * G REAL (KIND=R4KIND), PARAMETER :: PRT = 1.0 REAL (KIND=R4KIND), PARAMETER :: GAM1 = 0.2222222222222222222 !--------------------------------------------------------------------- !DCR: Nakanishi and Niino, 2009 REAL (KIND=R4KIND), PARAMETER :: A1 = 1.18 REAL (KIND=R4KIND), PARAMETER :: A2 = 0.665 REAL (KIND=R4KIND), PARAMETER :: B1 = 24.0 REAL (KIND=R4KIND), PARAMETER :: B2 = 15.0 REAL (KIND=R4KIND), PARAMETER :: C1 = 0.137 !------------------------------------------------------------------------- !J2001 REAL (KIND=R4KIND), PARAMETER :: A1 = 0.659888514560862645 !J2001 REAL (KIND=R4KIND), PARAMETER :: A2 = 0.6574209922667784586 !J2001 REAL (KIND=R4KIND), PARAMETER :: B1 = 11.87799326209552761 !J2001 REAL (KIND=R4KIND), PARAMETER :: B2 = 7.226971804046074028 !J2001 REAL (KIND=R4KIND), PARAMETER :: 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) !--------------------------------------------------------- REAL (KIND=R4KIND), PARAMETER :: RB1 = 1. / B1 !------------------------------------------- ! COEFFICIENTS OF THE TERMS IN THE NUMERATOR !------------------------------------------- REAL (KIND=R4KIND), PARAMETER :: ANMM = -3.*A1*A2*(3.*A2+3.*B2*C1+18.*A1*C1-B2)*BTG REAL (KIND=R4KIND), PARAMETER :: ANMH = -9.*A1*A2*A2*BTG*BTG REAL (KIND=R4KIND), PARAMETER :: BNMM = A1*(1.-3.*C1) REAL (KIND=R4KIND), PARAMETER :: BNMH = -A2*BTG !--------------------------------------------- ! COEFFICIENTS OF THE TERMS IN THE DENOMINATOR !--------------------------------------------- REAL (KIND=R4KIND), PARAMETER :: ADNM = 18.*A1*A1*A2*(B2-3.*A2)*BTG REAL (KIND=R4KIND), PARAMETER :: ADNH = 9.*A1*A2*A2*(12.*A1+3.*B2)*BTG*BTG REAL (KIND=R4KIND), PARAMETER :: BDNM = 6.*A1*A1 REAL (KIND=R4KIND), PARAMETER :: BDNH = 3.*A2*(7.*A1+B2)*BTG !----------------------------------------- ! COEFFICIENTS OF THE EQUILIBRIUM EQUATION !----------------------------------------- REAL (KIND=R4KIND), 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 (KIND=R4KIND), PARAMETER :: AEQH = 9.*A1*A2*A2*B1*BTG*BTG+9.*A1*A2*A2*(12.*A1+3.*B2) & & *BTG*BTG ! REAL (KIND=R4KIND), PARAMETER :: BEQM =-A1*B1*(1.-3.*C1)+6.*A1*A1 REAL (KIND=R4KIND), PARAMETER :: BEQH = A2*B1*BTG+3.*A2*(7.*A1+B2)*BTG !-------------------------- ! FORBIDDEN TURBULENCE AREA !-------------------------- REAL (KIND=R4KIND), PARAMETER :: REQU = -AEQH/AEQM*1.02 REAL (KIND=R4KIND), PARAMETER :: EPSGH =1.E-9 ! REAL (KIND=R4KIND), DIMENSION(LM) , INTENT(INOUT) ::& & Q2 ! REAL (KIND=R4KIND), DIMENSION(LM1) , INTENT(IN) ::& & GM , GH ! REAL (KIND=R4KIND), DIMENSION(LM1) , INTENT(INOUT) ::& & EL !------------------------ ! IMPLICIT NONE VARIABLES !------------------------ INTEGER(KIND=I4KIND) , INTENT(IN) ::& & LMHK ! REAL (KIND=R4KIND) , INTENT(IN) ::& & DTQ2 , USTAR ! INTEGER(KIND=I4KIND) ::& & LMHM , L ! REAL (KIND=R4KIND) ::& & GML , GHL , AEQU , BEQU , EQOL2 , ANUM , BNUM , ADEN , BDEN , & & CDEN , ARHS , BRHS , CRHS , DLOQ1 , ELOQ21 , ELOQ11 , ELOQ31 , ELOQ41 , & & ELOQ51 , RDEN1 , RHSP1 , ELOQ12 , ELOQ22 , ELOQ32 , ELOQ42 , ELOQ52 , RDEN2 , & & RHS2 , RHSP2 , RHST2 , ELOQ13 , ELOQN ! 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 END IF !------------------------ ! END OF TURBULENT BRANCH !------------------------ END IF !----------------------------------- ! 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