cGxke.for c c C.... File name .................... GXKE.F ..................... 040209 c

c
C#### JCL/MRM 20.12.12 check for GRND3 wall functions and don't
C####                  add GXTSKE source
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 '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/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#### SCM 04.02.09 Split commons for REALs and INTEGERs
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
c#### mrm 19.05.95 gxtske: provision for storing kt/kp & et/ep
      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 )
C#### mrm 25.07.02 Facility to modify k-e model constants
        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)
C#### mrm 25.07.02 Facility to modify k-e model constants
        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')
c#### mrm 19.05.95 gxtske: provision for storing kt/kp & et/ep
        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
c#### mrm 19.05.95 gxtske: provision for storing kt/kp & et/ep
            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)
C#### JCL/MRM 20.12.12 check for GRND3 wall functions and don't
C####                  add source
            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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/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