cGxke.for
c
c
C.... File name .................... GXKE.F ..................... 040209
c
c
c----------------------------------------------------------------------
c This file contains the following advanced-turbulence-option
c subroutines: GXREYT, GXREYN, GXLRDF, GXEYAP, GXFMU, GXFEPS,
c GXTSKE, GXRNGM,
c
C SUBROUTINE GXREYT is called from subroutine GXENUT and from
C subroutine GXKESO.
C
C
SUBROUTINE GXREYT(K1,K2,K3,K4)
C.... K1=REYT K2=KE K3=EP K4=VISL
include 'farray'
COMMON/RDATA/TINY,RGFIL1(84)
COMMON/IGE/IXF,IXL,IYF,IYL,IGFILL(21)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB='GXREYT'
IF(IXL.EQ.0) RETURN
CALL L0F4(K1,K2,K3,K4,I,I2M1,I3M1,I4M1,IADD,'GXREYT')
DO 1 IX=IXF,IXL
I=I+IADD
DO 1 IY=IYF,IYL
I=I+1
1 F(I)=AMAX1(F(I2M1+I)**2/((F(I3M1+I))*F(I4M1+I)+TINY),1.E-9)
NAMSUB='gxreyt'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXREYN is called from subroutine GXENUT and from
C subroutine GXKESO.
C
SUBROUTINE GXREYN(K1,K2,K3,K4)
C.... K1=REYN K2=KE K3=DIST K4=VISL
include 'farray'
COMMON/IGE/IXF,IXL,IYF,IYL,IGFILL(21)
COMMON/RDATA/TINY,RGFIL1(84)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB='GXREYN'
IF(IXL.EQ.0) RETURN
CALL L0F4(K1,K2,K3,K4,I,I2M1,I3M1,I4M1,IADD,'GXREYN')
DO 1 IX=IXF,IXL
I=I+IADD
DO 1 IY=IYF,IYL
I=I+1
1 F(I)=SQRT(ABS(F(I2M1+I))) * F(I3M1+I) / (F(I4M1+I)+TINY)
NAMSUB='gxreyn'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXLRDF is called from subroutine GXENUT and from
C subroutine GXKESO.
C
SUBROUTINE GXLRDF(JC,L0FMU,L0FONE,L0FTWO,L0REYT,L0REYN)
include 'farray'
COMMON/GENI/NXNY,IGFIL(59)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB='GXLRDF'
Compute Fmu
DO 2 J=1,NXNY
GREYN=F(L0REYN+J)
GREYT=F(L0REYT+J)
F(L0FMU+J)=1.0
IF(abs(GREYN).LT.800.0) F(L0FMU+J)=1.-EXP(-0.0165*GREYN)
F(L0FMU+J)=AMAX1(AMIN1(F(L0FMU+J)**2*(1.+20.5/GREYT),1.0),1.E-9)
2 CONTINUE
IF(JC.EQ.1) THEN
Compute Fone and Ftwo
DO 3 J=1,NXNY
F(L0FONE+J)=AMIN1(1.+(0.05/F(L0FMU+J))**3,1.E6)
F(L0FTWO+J)=1.0
GREYT=F(L0REYT+J)
IF(GREYT.LT.4.0) THEN
GRETSQ=AMIN1(GREYT*GREYT,100.)
F(L0FTWO+J)=AMAX1(1.-EXP(-GRETSQ),1.E-10)
ENDIF
3 CONTINUE
ENDIF
NAMSUB='gxlrdf'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXEYAP is called from group 13 OF GREX3 by setting
C the coefficient to FIXFLU and the value to GRND5 in the COVAL
C statement, and is entered only when the patch name begins with
C the character 'EYAP' and the variable LTLS is solved to
C calculate the distance to the nearest wall.
C
SUBROUTINE GXEYAP
include 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
INTEGER VAL,CO,WALDIS,PATGEO
COMMON/RDATA/TINY,RDATS(84)
COMMON/IDATA/NX,NY,IFILL(118)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXEYAP'
CALL SUB3(L0VAL,L0F(VAL),L0KE,L0F(KE),L0EP,L0F(EP))
L0WDIS=L0F(LBNAME('WDIS'))
GYAPC=AK/CMUCD**0.75
DO 1 I=1,NX*NY
IF(F(L0WDIS+I).LT.TINY) THEN
F(L0VAL+I)=0.0
ELSE
AKE=AMAX1(1.E-10,F(L0KE+I))
AEP=AMAX1(1.E-10,F(L0EP+I))
GLENRT=(AKE**1.5/AEP)/(GYAPC*F(L0WDIS+I))
GYAPS =0.83*(GLENRT-1.0)*GLENRT*GLENRT*AEP**2/AKE
F(L0VAL+I)=AMAX1(0.0,GYAPS)
ENDIF
1 CONTINUE
NAMSUB = 'gxeyap'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXFMU computes the ENUT damping function required
C by the two-layer k-eps model. The subroutine is called from
C subroutines GXENUT and GXKESO.
C
SUBROUTINE GXFMU(L0FMU,L0REYN)
include 'farray'
COMMON/GENI/NXNY,IGFIL(59)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB='GXFMU '
DO 2 J=1,NXNY
GREYN=F(L0REYN+J)
F(L0FMU+J)=1.0
IF(ABS(GREYN).LT.350.0)
1 F(L0FMU+J)=AMAX1(AMIN1(1.-EXP(-0.0198*GREYN),1.0),1.E-9)
2 CONTINUE
NAMSUB='gxfmu '
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXFEPS computes the EP damping function required
C by the two-layer k-eps model. The subroutine is called from
C subroutine GXKESO.
C
SUBROUTINE GXFEPS(L0FTWO,L0REYN)
include 'farray'
COMMON/GENI/NXNY,IGFIL(59)
COMMON/RDATA/TINY,RDATS(84)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB='GXFEPS'
DO 2 J=1,NXNY
GREYN=F(L0REYN+J)
F(L0FTWO+J)=1.0
IF(ABS(GREYN).LT.350.0)
1 F(L0FTWO+J)=AMIN1(1.E6,1.+5.3/(GREYN+TINY))
2 CONTINUE
NAMSUB='gxfeps'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXTSKE is called from group 13 OF GREX3 for variables:
C EP, KP, ET, and KT, in that order, when the coefficient & value
C are set to GRND5 in the COVAL statement, and the name of a patch,
C of PHASEM type, begins with the characters 'TSKE'.
C
SUBROUTINE GXTSKE(VIST,LEN1,FIXFLU)
include 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON/RDA4/PRT(150)
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
INTEGER VAL,CO,WALDIS,PATGEO,VIST
COMMON/RDATA/TINY,RGFIL(84)
COMMON/IDATA/NX,NY,IFILL(118)
C COMMON/TSKEM/ GKTDKP,LBKP,LBKT,LBET,LBVOSQ,LBOMEG
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
1 /TSKEMR/ GKTDKP
COMMON/NAMFN/NAMFUN,NAMSUB
COMMON/SORSM/SORSUM
CHARACTER*6 NAMFUN,NAMSUB
INTEGER COGRN,VALGRN
LOGICAL ZRGRN3
SAVE NXNY,GC3E,GCT1,GCT2,GCT3,L0PRAT,L0TRAT,L0EPOL,L0KE,L0EP,
1 L0PROD,LBKTKP,LBETEP
C
NAMSUB = 'GXTSKE'
IF(IGR.EQ.1 .AND. ISC.EQ.3) THEN
CALL SUB4R( GC3E,0.21, GCT1,0.29, GCT2,1.28, GCT3,1.66 )
CALL GETSPD('KECONST','GC3E',1,GC3E,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','GCT1',1,GCT1,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','GCT2',1,GCT2,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('KECONST','GCT3',1,GCT3,IDUM,.FALSE.,' ',IERR)
CALL WRITBL
CALL WRIT40('2-scale K-E turbulence model constants ')
CALL WRIT2R('GC3E ',GC3E ,'GCT1 ',GCT1)
CALL WRIT2R('GCT2 ',GCT2 ,'GCT3 ',GCT3)
CALL WRITBL
NXNY=NX*NY
CALL SUB3( LBKP,LBNAME('KP '), LBKT,LBNAME('KT '),
1 LBET,LBNAME('ET ') )
CALL GXMAKE(L0PRAT,NXNY,'PRAT')
CALL GXMAKE(L0TRAT,NXNY,'TRAT')
CALL GXMAKE(L0EPOL,NXNY,'EPOL')
CALL GXMAKE(L0PROD,NXNY,'PROD')
CALL SUB2(LBKTKP,LBNAME('KTKP'),LBETEP,LBNAME('ETEP'))
ELSEIF(IGR.EQ.13) THEN
CALL SUB2(COGRN,ISC-1,VALGRN,ISC-12)
CALL SUB3(L0KE,L0F(KE),L0KT,L0F(LBKT),L0ET,L0F(LBET))
C
IF(INDVAR.EQ.EP) THEN
c ep, energy transfer rate from production to dissipation range
c..... this section is entered first at each slab
IF(COGRN.EQ.5) THEN
L0CO=L0F(CO)
L0KP=L0F(LBKP)
L0EP=L0F(EP)
L0KT=L0F(LBKT)
L0KE=L0F(KE)
L0GEN1=L0F(GEN1)
L0VIST=L0F(VIST)
L0M1=L0F(LM1)
FACTOR=1.3333*C2E
DO 3 I=1,NXNY
F(L0KE+I)=F(L0KP+I)+F(L0KT+I)
c.... create the rate factor for "production" quantities
c store the start of sweep ep
c store the production rate
F(L0PRAT+I)=F(L0EP+I)/(F(L0KP+I)+TINY)
F(L0EPOL+I)=F(L0EP+I)
F(L0PROD+I)=F(L0GEN1+I)*F(L0VIST+I)
c.... create the coefficient
F(L0CO+I)=F(L0PRAT+I)*FACTOR
3 CONTINUE
IF(LBKTKP.NE.0) CALL FN15(LBKTKP,LBKT,LBKP,0.0,1.0)
IF(LBETEP.NE.0) CALL FN15(LBETEP,LBET,EP,0.0,1.0)
ELSEIF(VALGRN.EQ.5) THEN
L0VAL=L0F(VAL)
L0EP=L0F(EP)
L0KE=L0F(KE)
L0SU=L0F(LSU)
L0M1=L0F(LM1)
DO 4 I=1,NXNY
F(L0VAL+I)=F(L0EPOL+I)*0.25
SU=F(L0PROD+I)*F(L0PRAT+I)*F(L0M1+I)
1 * (C1E+GC3E*F(L0PROD+I)/(F(L0EPOL+I)+TINY))
F(L0SU+I)=F(L0SU+I)+SU
SORSUM=SORSUM+SU
4 CONTINUE
ENDIF
C
ELSEIF(INDVAR.EQ.LBKP) THEN
c kp, turbulence energy in production range
IF(COGRN.EQ.5) THEN
L0CO=L0F(CO)
DO 1 I=1,NXNY
F(L0CO+I)=F(L0PRAT+I)*1.5
1 CONTINUE
ELSEIF(VALGRN.EQ.5) THEN
L0VAL=L0F(VAL)
L0SU=L0F(LSU)
L0KP=L0F(LBKP)
L0M1=L0F(LM1)
DO IX=1,NX
DO IY=1,NY
I=(IX-1)*NY+IY
C.... Wall-type patch is present in the cell, so loop over wall-patches
C to check for GRND3 coefficient for KP
IF(ZRGRN3(IX,IY,IZSTEP,I)) THEN
F(L0VAL+I)=F(L0KP+I)
SU = 0.0
ELSE
F(L0VAL+I)=0.3333*F(L0KP+I)
SU= F(L0PROD+I)*F(L0M1+I)
ENDIF
F(L0SU+I)=F(L0SU+I)+SU
SORSUM=SORSUM+SU
ENDDO
ENDDO
ENDIF
C
ELSEIF(INDVAR.EQ.LBET) THEN
c et, dissipation rate of kt
L0ET=L0F(LBET)
L0KT=L0F(LBKT)
IF(COGRN.EQ.5) THEN
L0CO=L0F(CO)
DO 7 I=1,NXNY
F(L0TRAT+I)=F(L0ET+I)/(F(L0KT+I)+TINY)
F(L0CO+I) =F(L0TRAT+I)*GCT3*1.3333
7 CONTINUE
ELSEIF(VALGRN.EQ.5) THEN
L0VAL=L0F(VAL)
L0SU=L0F(LSU)
L0M1=L0F(LM1)
DO 6 I=1,NXNY
GEP=F(L0EPOL+I)
F(L0VAL+I)=0.25*F(L0ET+I)
SU=F(L0M1+I)*(GCT2*F(L0ET+I) +
1 GCT1*GEP )*GEP/(F(L0KT+I)+TINY)
F(L0SU+I)=F(L0SU+I)+SU
SORSUM=SORSUM+SU
6 CONTINUE
ENDIF
ELSEIF(INDVAR.EQ.LBKT) THEN
c kt, turbulence energy in dissipation range
IF(COGRN.EQ.5) THEN
CALL FN2(CO,-L0TRAT,0.0,1.5)
ELSEIF(VALGRN.EQ.5) THEN
L0VAL=L0F(VAL)
L0SU=L0F(LSU)
L0M1=L0F(LM1)
L0KT=L0F(LBKT)
DO 5 I=1,NXNY
F(L0VAL+I) = 0.3333*F(L0KT+I)
SU= F(L0EPOL+I)*F(L0M1+I)
F(L0SU+I)=F(L0SU+I)+SU
SORSUM=SORSUM+SU
5 CONTINUE
ENDIF
ENDIF
ENDIF
NAMSUB = 'gxtske'
END
c file-name gxrng.for 220394
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXRNGM is called from group 13 OF GREX3 by setting
C the coefficient & value to GRND4 in the COVAL statement, and
C is entered only when the patch name begins with the
C character 'RNGM'.
C
C.... The arguments VIST and LEN1 are the integer names used in the
C SATELLITE, indicating which whole-field store will be used for
C the turbulent kinematic viscosity and length scale of phase 1
C respectively.
C
C.... The library cases T300-T306 make use of it.
C
C
SUBROUTINE GXRNGM(VIST,LEN1,FIXFLU)
include 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
INTEGER VAL,CO,WALDIS,PATGEO,VIST
COMMON/IDATA/NX,NY,IFILL(118)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
INTEGER COGRN,VALGRN
SAVE ETA0,BETA,NXNY,JALF,JETA,L0ALF,L0ETA,L0KE,L0EP
C
NAMSUB = 'GXRNGM'
IF(IGR.EQ.1 .AND. ISC.EQ.3) THEN
CALL SUB2R( ETA0,4.38, BETA,0.012 )
NXNY= NX*NY
CALL SUB2( JALF,LBNAME('ALF '), JETA,LBNAME('ETA ') )
IF(JALF.EQ.0) CALL GXMAKE(L0ALF,NXNY,'ALF ')
IF(JETA.EQ.0) CALL GXMAKE(L0ETA,NXNY,'ETA ')
ELSE
CALL SUB2(COGRN,ISC-1,VALGRN,ISC-12)
IF(COGRN.EQ.4) THEN
C.... Coefficient part of the linearized rate-of-strain term R ------CO
C.... If alf > 0, CO = ALF*EP/KE otherwise CO = FIXFLU
C
C where ALF = CMUCD*ETA**3*(1-ETA/ETA0)/(1+BETA*ETA**3) and
C ETA = SQRT(GEN1)*KE/EP
C
CALL SUB3(L0CO,L0F(CO),L0KE,L0F(KE),L0EP,L0F(EP))
L0GEN1=L0F(GEN1)
IF(JALF.NE.0) L0ALF=L0F(JALF)
IF(JETA.NE.0) L0ETA=L0F(JETA)
DO 1 I=1,NXNY
rateco=f(l0ep+i)/(f(l0ke+i)+1.e-20) +1.e-10
F(L0ETA+I)=amin1(1.e3, SQRT(abs(F(L0GEN1+I)))/rateco)
ETA3=F(L0ETA+I)**3
F(L0ALF+I)=CMUCD*ETA3*(ETA0-F(L0ETA+I))/
1 (ETA0*(1.0+BETA*ETA3))
IF(F(L0ALF+I).LT.0) THEN
F(L0CO+I)=FIXFLU
ELSE
F(L0CO+I)=F(L0ALF+I)*rateco
ENDIF
1 CONTINUE
ELSEIF(VALGRN.EQ.4) THEN
C.... Value part of the linearized source ---------------------------VAL
C.... If alf > 0, VAL = 0.0 otherwise VAL = -ALF*EP**2/K
C
L0VAL=L0F(VAL)
FAC=1.0/FIXFLU
DO 2 I=1,NXNY
IF(F(L0ALF+I).LT.0) THEN
F(L0VAL+I)=-F(L0ALF+I)*F(L0EP+I)*F(L0EP+I)*FAC/
1 (F(L0KE+I)+1.E-20)
ELSE
F(L0VAL+I)=0.0
ENDIF
2 CONTINUE
ENDIF
ENDIF
NAMSUB = 'gxrngm'
END
c