c

C
      REAL FUNCTION KINVIST(I)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/mfmdat'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
     1            ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
     1      /F0/IF01(29),L0XYDX,L0XYDY,IF02(4),L0XYXG,IF03,L0XYYG,IF04,
     1          L0XYDXG,L0XYDYG,IF05(124),L0XYDZ,L0XYDZG,IF06(137)
     1      /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
     1              L0UD2,L0UD3,L0UD4
     1      /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
     1      /LRNTM3/L0UTAU,NUMWAL,L0WALL,L0DSKN,IVPRST
     1      /MMKKE/ LBFOMG,L0FOMG
     1      /DOMSIZ/SIZMAX
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
     1      /MCRMXI/LBRATE,L0APXY,L0SUXY,L0APCL,L0SUCL,L0RATE
C
C.... KINVIST is turbulent contribution to total kinematic viscosity
      IF(IGRND.EQ.-1) THEN           ! use the PIL variable ENUT
        KINVIST=PRPRTY
      ELSEIF(IGRND.EQ.1) THEN        ! linear function of length
        KINVIST= ENUTA + ENUTB*F(L0LEN+I)
      ELSEIF(IGRND.EQ.2) THEN        ! Prandtl mixing-length formula
        KINVIST= F(L0LEN+I)**2 * SQRT(ABS(F(L0GEN+I)))
      ELSEIF(IGRND.EQ.3) THEN        ! Prandtl-Kolmogorov formula
        KINVIST= CMU*F(L0LEN+I)*SQRT(ABS(F(L0KE+I)))
      ELSEIF(IGRND.EQ.4) THEN        ! two-fluid model of turbulence
        ABSDV = VSLPCL(I,XCYCZ,SLUNX1)
        KINVIST= ENUTA*F(L0LEN+I)*ABSDV*F(L0R1+I)*F(L0R2+I)
      ELSEIF(IGRND.EQ.5) THEN        ! Harlow-Nakayama formula
        KINVIST= CMUCD*F(L0KE+I)**2/AMAX1(F(L0EP+I),TINY)
      ELSEIF(IGRND.EQ.6) THEN        ! Saffman-Spalding formula
        KINVIST= CMUCD*F(L0KE+I)/SQRT(AMAX1(F(L0VOSQ+I),TINY))
      ELSEIF(IGRND.EQ.7) THEN        ! Kolmogorov formula
        KINVIST= F(L0KE+I)/AMAX1(F(L0OMEG+I),TINY)
      ELSEIF(IGRND.EQ.8) THEN        ! LVEL model (SKNFG2 returns the
        VISLAM= F(L0VISL+I)          ! skin-friction factor at the wall)
        REYN  = SQRT(ABSVSQ(I,XCYCZ,SLUNX1))*F(L0LEN+I)/(VISLAM+TINY)
        IF(REYN.LT.1.0) THEN         ! very close to the wall there
          KINVIST= 0.0               ! is no turbulent contribution
        ELSE                         ! formula from Spalding's wall law
          TAUP= SKNFG2(REYN,5+ISKINA,ISKINB,.FALSE.,REYN,EEFF,SKINT)
          AKUP= AK/SQRT(TAUP)        ! K times uplus
          SUM = ( (0.1666667*AKUP + 0.5) * AKUP + 1.0) * AKUP + 1.0
          KINVIST = VISLAM * ( EXP(AKUP) - SUM ) * AK / EWAL
        ENDIF
      ELSEIF(IGRND.EQ.9) THEN
        KINVIST= ENUTA               ! uniform value
      ELSEIF(IGRND.EQ.10) THEN
        KINVIST= VISCON*F(L0MIXL+I)**2*F(L0RATE+I) ! Formula used in some MFM
      ENDIF                                        ! cases
C.... Modifications for low Reynolds numbers:
      IF(IENUTA.EQ.3 .OR.IENUTA.EQ.4) THEN  ! Lam-Bremhorst K-E model
        FMU = 1.0
        REYT= AMAX1( F(L0KE+I)**2/(F(L0EP+I)*F(L0VISL+I)+TINY) , 1.E-9)
        REYN= SQRT(ABS(F(L0KE+I)))*F(L0WDIS+I)/(F(L0VISL+I)+TINY)
        IF(ABS(REYN).LT.800.) FMU= 1.0 - EXP(-0.0165*REYN)
        FMU= AMAX1(AMIN1(FMU**2*(1.+20.5/REYT), 1.0), 1.E-9)
        KINVIST= KINVIST*FMU
        F(L0REYT+I)= REYT
        F(L0REYN+I)= REYN
        F(L0FMU +I)= FMU
      ELSEIF(IENUTA.EQ.6) THEN  ! a Simple low-Re modification of ENUT
C     VisT = VisT_nom * amin1(1.0 , (A * VisT_nom / VisL ) ** B)
C     where:   VisT     = turbulent kinematic viscosity to be used,
C              VisT_nom = nominal (high-Reynolds-number) value of VisT
C              VisL     = laminar viscosity
C     Suitable values for A and B may be 0.2 and 4.0
        VISLAM= F(L0VISL+I)
        FACTOR= ENUTA*KINVIST
        IF(FACTOR.LT.VISLAM) KINVIST= KINVIST*(FACTOR/VISLAM)**ENUTB
      ELSEIF(IENUTA.EQ.8) THEN  ! The 2-layer model
        FMU = 1.0
        REYN= SQRT(ABS(F(L0KE+I)))*F(L0WDIS+I)/(F(L0VISL+I)+TINY)
        IF(ABS(REYN).LT.350.) THEN
          FMU= AMAX1(AMIN1(1.0 - EXP(-0.0198*REYN), 1.0), 1.E-9)
          KINVIST= KINVIST*FMU
        ENDIF
        F(L0REYN+I)= REYN
        F(L0FMU +I)= FMU
      ELSEIF(IENUTA.EQ.9) THEN  ! Another simple low-Re modification
C     when Const*VisT < VisL, VisT is set equal to zero.
C     A suitable value for ENUTA may be 0.05
        IF(ENUTA*KINVIST.LT.F(L0VISL+I)) KINVIST= 0.0
      ELSEIF(IENUTA.EQ.11) THEN ! for the K-Omega model
        REYT= KINVIST/(F(L0VISL+I)+TINY)
        FMU = (0.15 + REYT)/(6.0 + REYT)
        KINVIST= KINVIST*FMU
        F(L0REYT+I)= REYT
        F(L0FMU +I)= FMU
      ELSEIF(IENUTA.EQ.12) THEN
C.... MMK high-Reynolds-number K-E turbulence model
        FOMEG = F(L0FOMG+I)
        IF(FOMEG.LT.1.0) KINVIST= KINVIST*FOMEG
      ENDIF
c..................................... end of low-Re modifications
      IF(.NOT.ONEPHS) THEN    ! for two-phase flow
C.... Bubble-induced contribution to turbulent kinematic viscosity:
C       Enut= Enut + ENUTC * R1 * R2 * Vslp * (possibly) len
        IF(ENUTC.GT.0.0) THEN
          VSLP  = AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1))
          VISADD = ENUTC * F(L0R2+I) * (1.0 - F(L0R2+I)) * VSLP
          IF(L0LEN.NE.0) THEN
            VISADD = VISADD * F(L0LEN+I)
          ENDIF
          KINVIST = KINVIST + VISADD
        ENDIF
      ENDIF
      END
C***********************************************************************
      SUBROUTINE SLBVST(IPILOPT,dbgloc)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/VMSCMN/FL1CON
     1      /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
     1              L0UD2,L0UD3,L0UD4
     1      /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
     1      /MMKKE/ LBFOMG,L0FOMG
     1      /VELCMN/L0UVW(6),L0UVW2(6) /FLPCMN/IFILP(30)
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
     1      /TSKEMR/ GKTDKP
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,NWHOLE, IGSH,
     1           IGF3(19),IPRL,
     1           IBTAU,IGF4(16),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
     1           IPRPS,IRADX,IRADY,IRADZ,IVFOL
     1      /F0/ IF01(29),L0XYDX,L0XYDY,IF02(3),L0XYRV,L0XYXG,IF03,
     1           L0XYYG,IF04,L0XYDXG,L0XYDYG,IF05(68),KZXCY,IF06(37),
     1           L0AHZ,IF07(17),L0XYDZ,L0XYDZG,IF08(137)
      COMMON/NAMFN/NAMFUN,NAMSUB
      REAL KINVIST
      LOGICAL DBGLOC,SLD,FL1CON,GRN
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB= 'SLBVST'
      if(flag.or.dbgloc) call banner(1,'SLBVST',280311 )
        IGR= 9
        ISC= IPROP
C.... Call GROUND for user-set property:
      IF(IPILOPT.EQ.0) THEN
        IF(USEGRD) THEN
          CALL GROUND
        ENDIF
        GO TO 800
      ENDIF
      IF(IPILOPT.EQ.-1) GO TO 700
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
      IF(IGRND.EQ.4.OR.IGRND.EQ.8.OR.(.NOT.ONEPHS.AND.ENUTC.GT.0.0))THEN
        IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
        IF(SOLVE(3)) L0UVW(3)=L0F(3)
        IF(SOLVE(5)) L0UVW(1)=L0F(5)
        IF(SOLVE(7)) L0UVW(5)=L0F(7)
        IF(.NOT.ONEPHS) THEN
          IF(SOLVE(4)) L0UVW2(3)=L0F(4)
          IF(SOLVE(6)) L0UVW2(1)=L0F(6)
          IF(SOLVE(8)) L0UVW2(5)=L0F(8)
        ENDIF
      ENDIF
      L0VISL  = L0F(VISL)
      IF(IGRND.EQ.3) THEN
        IF(STORE(12)) L0KE=L0F(12)
      ELSEIF(IGRND.EQ.4) THEN
        IF(STORE(9))  L0R1=L0F(9)
        IF(SOLVE(10)) L0R2=L0F(10)
      ELSEIF(IGRND.EQ.5) THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(STORE(13)) L0EP=L0F(13)
      ELSEIF(IGRND.EQ.6)THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(LBVOSQ.GT.0) L0VOSQ=L0F(LBVOSQ)
      ELSEIF(IGRND.EQ.7) THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(LBOMEG.GT.0) L0OMEG=L0F(LBOMEG)
      ELSEIF(IGRND.EQ.8) THEN
        L0VISL=L0F(VISL)
        IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
      ENDIF
      IF(USEWDS) L0WDIS=L0F(LBWDIS)
      IF(GENK) L0GEN=L0F(GEN1)
      IF(LEN1.GT.0 .OR. GRN(EL1)) L0LEN=L0F(LEN1)
      IF(IENUTA.EQ.3 .OR.IENUTA.EQ.4) THEN
        IF(LBREYT.NE.0) L0REYT= L0F(LBREYT)
        IF(LBREYN.NE.0) L0REYN= L0F(LBREYN)
        IF(LBFMU.NE.0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA.EQ.8) THEN
        IF(LBREYN.NE.0) L0REYN= L0F(LBREYN)
        IF(LBFMU.NE.0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA.EQ.11) THEN
        IF(LBREYT.NE.0) L0REYT= L0F(LBREYT)
        IF(LBFMU.NE.0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA.EQ.12) THEN
        IF(LBFOMG.NE.0)  L0FOMG = L0F(LBFOMG)
      ENDIF
C-----------------------------------------------------------------------
C.... Loop over slab to get and set cell properties:
 700  IGRND=IPILOPT
      IF(IPRPS.EQ.0) THEN
C.... One material only
        DO 60 I= 1,NXNYST
  60    F(KPROP+I)= KINVIST(I)
      ELSE
C.... exclude solids
        DO 70 I= 1,NXNYST
          IF(SLD(I)) THEN
            F(KPROP+I)= TINY
          ELSE
            F(KPROP+I)= KINVIST(I)
          ENDIF
  70    CONTINUE
      ENDIF
C----------------------------------------------------------------------
C.... Corrections, debug print-out, and other property adjustments:
      IF(IPILOPT.EQ.8 .AND. LBEPKE.NE.0) THEN
        L0EPKE= L0F(LBEPKE)
        DO 80 I= 1,NXNYST
          IF(SLD(I)) THEN
            F(L0EPKE+I)= 0.0
          ELSE
            RWDIS= 1./(F(L0WDIS+I)+TINY)
            F(L0EPKE+I)= F(KPROP+I)*RWDIS*RWDIS
          ENDIF
  80    CONTINUE
      ENDIF
C.... Call GREX to correct a property set above:
 800  IF(USEGRX) CALL GREX3
C.... Call ALTPRP for an alternative property setting
      IF(USEALT) CALL ALTPRP
C.... Call GROUND for the user to correct a property set above:
      IF(USEGRD) THEN
        IF(IPILOPT.GT.0) CALL GROUND
      ENDIF
      if(flag.or.dbgloc) call banner(2,'SLBVST',280311 )
      NAMSUB= 'SLBVST'
      END
c-------------------------
      SUBROUTINE INIVST
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/qdata'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/prpcmn'
      COMMON/TURB/CMU,CD,CMUCD,C1E,C2E,AK,EWAL
     1      /SPEPRC/MFM
     1      /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
     1              L0UD2,L0UD3,L0UD4/GENI/IGN1(49),ITEM1,ITEM2,IGN2(9)
     1      /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
     1      /LRNTM3/L0UTAU,NUMWAL,L0WALL,L0DSKN,IVPRST
     1      /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
      COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
     1      /TSKEMR/ GKTDKP
     1      /VELCMN/L0UVW(6),L0UVW2(6) /FLPCMN/IFILP(30)
      LOGICAL GETREALS
      DIMENSION NAMPRP(30)
      EQUIVALENCE (NAMPRP(1),DENST1)
C.... Preliminaries:
      if(flag) call banner(1,'INIVST',280311 )
C.... Turbulence-model constants:
      IF(.NOT.RSTM) THEN
        IF(CMU.LE.0.0)  CMU=0.5478
        IF(CD.LE.0.0)   CD= 0.1643
        IF(C1E.LE.0.0)  C1E=1.44
        IF(C2E.LE.0.0)  C2E=1.92
        IF(AK.LE.0.0)   AK= 0.41
        IF(EWAL.LE.0.0) EWAL=8.6
c
        IF(GETREALS('TURCONST')) THEN
          CMU=REALS(1) ; CD=REALS(2) ; C1E=REALS(3)
          C2E=REALS(4) ; AK=REALS(5) ; EWAL=REALS(6)
        ELSE
          CALL GETSPD('KECONST','CMU',1,CMU,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CD',1,CD,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','AK' ,1,AK,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','EWAL' ,1,EWAL,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
        ENDIF
        CMUCD=CMU*CD
        RSCMCD=CMUCD**(-0.5) ; CDDAK2=2.0*CD/AK ; TAUDKE=SQRT(CMUCD)
        RTTDKE=SQRT(TAUDKE) ; AKC=AK*RTTDKE ; EWC=EWAL*RTTDKE
        IF(IENUTA.EQ.1) THEN
C.... RNG-derived K-E model:
          CMUCD=0.0845 ; C1E=1.42 ; C2E=1.68
          CD = CMUCD**0.75 ; CMU= CMUCD/CD
          CALL GETSPD('KECONST','CMU',1,CMU,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CD',1,CD,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CMUCD=CMU*CD
          CALL GXRNGM(0,0,FIXFLU)
        ELSEIF(IENUTA.EQ.2 .OR. IENUTA.EQ.4) THEN
C.... Chen-Kim high & low-Re K-E models
          C1E=1.15 ; C2E=1.9
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
        ELSEIF(IENUTA.EQ.7) THEN
C.... Two-scale K-E model:
          C1E=1.24 ; C2E=1.84 ; GC3E=0.21
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','GC3E',1,GC3E,IDUM,.FALSE.,'   ',IERR)
          GKTDKP= AK*AK/(PRT(13)*SQRT(CMUCD)*(C2E-GC3E-C1E)) - 1.0
          CALL GXTSKE(0,0,FIXFLU)
        ELSEIF(IENUTA.EQ.10 .OR. IENUTA.EQ.11) THEN
C.... Wilcox-Kolmogorov high & low-Re K-W models
          C1E=5./9. ; C2E=0.075
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          KELIN= 0
        ENDIF
        IF(.NOT.NULLPR) THEN
          CALL WRITBL
          CALL WRIT40('>> Current turbulence model constants <<')
          CALL WRIT40('They may be changed by inserting in Q1  ')
          CALL WRIT40('SPEDAT(KECONST,name of constant,R,value)')
          CALL WRITBL
          CALL WRIT3R('CMU     ',CMU     ,
     1                'CD      ',CD,'CMUCD   ',CMUCD)
          CALL WRIT2R('C1E     ',C1E     ,'C2E     ',C2E)
          CALL WRIT2R('AK      ',AK      ,'EWAL    ',EWAL)
          CALL WRITBL
        ENDIF
      ENDIF
      RSCMCD=CMUCD**(-0.5) ; CDDAK2=2.0*CD/AK ; TAUDKE=SQRT(CMUCD)
      RTTDKE=SQRT(TAUDKE) ; AKC=AK*RTTDKE ; EWC=EWAL*RTTDKE
C.... Memory allocation for use in turbulence modelling:
      IF(ENUT.LT.0.0) THEN
        CALL ADDL1C(L0UTAU,IENUTA.EQ.5,'l0utau')
        IF(IENUTA.EQ.3 .OR. IENUTA.EQ.4 .OR. IENUTA.EQ.8 .OR.
     1     IENUTA.EQ.11) THEN
          NXNY=NX*NY
          IF(LBFONE.EQ.0) CALL GXMAKE(L0FONE,NXNY,'FONE')
          IF(LBFTWO.EQ.0) CALL GXMAKE(L0FTWO,NXNY,'FTWO')
          IF(LBFMU.EQ.0)  CALL GXMAKE(L0FMU,NXNY,'FMU ')
          IF(LBREYN.EQ.0.AND.IENUTA.NE.11) CALL GXMAKE(L0REYN,NXNY,'RE')
          IF(LBREYT.EQ.0.AND.IENUTA.NE.8) CALL GXMAKE(L0REYT,NXNY,'RET')
        ENDIF
      ENDIF
      if(flag) call banner(2,'INIVST',280311 )
      END
c