c

C
      REAL FUNCTION KINVIST(I)
      INCLUDE 'farray'
      INCLUDE 'mfmdat'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE '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==-1) THEN           ! use the PIL variable ENUT
        KINVIST=PRPRTY
      ELSEIF(IGRND==1) THEN        ! linear function of length
        KINVIST= ENUTA + ENUTB*F(L0LEN+I)
      ELSEIF(IGRND==2) THEN        ! Prandtl mixing-length formula
        KINVIST= F(L0LEN+I)**2 * SQRT(ABS(F(L0GEN+I)))
      ELSEIF(IGRND==3) THEN        ! Prandtl-Kolmogorov formula
        KINVIST= CMU*F(L0LEN+I)*SQRT(ABS(F(L0KE+I)))
      ELSEIF(IGRND==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==5) THEN        ! Harlow-Nakayama formula
        KINVIST= CMUCD*F(L0KE+I)**2/AMAX1(F(L0EP+I),TINY)
      ELSEIF(IGRND==6) THEN        ! Saffman-Spalding formula
        KINVIST= CMUCD*F(L0KE+I)/SQRT(AMAX1(F(L0VOSQ+I),TINY))
      ELSEIF(IGRND==7) THEN        ! Kolmogorov formula
        KINVIST= F(L0KE+I)/AMAX1(F(L0OMEG+I),TINY)
      ELSEIF(IGRND==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(LBREYN>0) THEN
          F(L0REYN+I)=REYN
        ENDIF
        IF(REYN<1.0) THEN         ! very close to the wall there
          KINVIST= 0.0               ! is no turbulent contribution
        ELSE                         ! formula from Spalding's wall law
c          TAUP= SKNFG2(REYN,5+ISKINA,ISKINB,.FALSE.,REYN,EEFF,SKINT)
c          AKUP= AK/SQRT(TAUP)        ! K times uplus
c          SUM = ( (0.1666667*AKUP + 0.5) * AKUP + 1.0) * AKUP + 1.0
c          KINVIST = VISLAM * ( EXP(AKUP) - SUM ) * AK / EWAL
c dbs 17.06.12 The above recourse to sknfg2 is not necessary; and indeed
c              it is incorrect; for  there is no clearly-connected
c              near-surface cell where the skin friction can be calculated
c              for circumstances sketched here.
c
c                            x <-- point for which effective viscosity
c                                  is to be calculated
c       ////////////////
c       ///// solid ////
c       ////////////////
c
c              The following call to enutpl   suffice:
c
          KINVIST=VISLAM*ENUTPL(REYN)
        ENDIF
      ELSEIF(IGRND==9) THEN
        KINVIST= ENUTA               ! uniform value
      ELSEIF(IGRND==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==3 .OR.IENUTA==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)<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==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(FACTOR0.0) THEN
          VSLP  = AMAX1(CFIPA, VSLPCL(I,XCYCZ,SLUNX1))
          VISADD = ENUTC * F(L0R2+I) * (1.0 - F(L0R2+I)) * VSLP
          IF(L0LEN/=0) THEN
            VISADD = VISADD * F(L0LEN+I)
          ENDIF
          KINVIST = KINVIST + VISADD
        ENDIF
      ENDIF
      END
C***********************************************************************
      SUBROUTINE SLBVST(IPILOPT,dbgloc)
      INCLUDE 'farray'
      INCLUDE 'mfmdat'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'satear'
      INCLUDE 'grdear'
      INCLUDE '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 /MCRMXI/LBRATE,L0APXY,L0SUXY,L0APCL,L0SUCL,L0RATE,L0FRCS,
     1               L0ABSC,L0AV,L0BV,L0CV
      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',131112 )
        IGR= 9
        ISC= IPROP
C.... Call GROUND for user-set property:
      IF(IPILOPT==0) THEN
        IF(USEGRD) THEN
          CALL GROUND
        ENDIF
        GO TO 800
      ENDIF
      IF(IPILOPT==-1) GO TO 700
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
      IF(IGRND==4.OR.IGRND==8.OR.(.NOT.ONEPHS.AND.ENUTC>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==3) THEN
        IF(STORE(12)) L0KE=L0F(12)
      ELSEIF(IGRND==4) THEN
        IF(STORE(9))  L0R1=L0F(9)
        IF(SOLVE(10)) L0R2=L0F(10)
      ELSEIF(IGRND==5) THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(STORE(13)) L0EP=L0F(13)
      ELSEIF(IGRND==6)THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(LBVOSQ>0) L0VOSQ=L0F(LBVOSQ)
      ELSEIF(IGRND==7) THEN
        IF(STORE(12)) L0KE=L0F(12)
        IF(LBOMEG>0) L0OMEG=L0F(LBOMEG)
      ELSEIF(IGRND==8) THEN
        L0VISL=L0F(VISL)
        IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
        IF(LBREYN/=0) L0REYN= L0F(LBREYN)
      ELSEIF(IGRND==10) THEN
        IF(LBRATE>0) L0RATE=L0F(LBRATE)
        IF(LBMIXL>0) L0MIXL=L0F(LBMIXL)
      ENDIF
      IF(USEWDS) L0WDIS=L0F(LBWDIS)
      IF(GENK) L0GEN=L0F(GEN1)
      IF(LEN1>0 .OR. GRN(EL1)) L0LEN=L0F(LEN1)
      IF(IENUTA==3 .OR.IENUTA==4) THEN
        IF(LBREYT/=0) L0REYT= L0F(LBREYT)
        IF(LBREYN/=0) L0REYN= L0F(LBREYN)
        IF(LBFMU/=0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA==8) THEN
        IF(LBREYN/=0) L0REYN= L0F(LBREYN)
        IF(LBFMU/=0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA==11) THEN
        IF(LBREYT/=0) L0REYT= L0F(LBREYT)
        IF(LBFMU/=0)  L0FMU = L0F(LBFMU)
      ELSEIF(IENUTA==12) THEN
        IF(LBFOMG/=0)  L0FOMG = L0F(LBFOMG)
      ENDIF
C-----------------------------------------------------------------------
C.... Loop over slab to get and set cell properties:
 700  IGRND=IPILOPT
      IF(IPRPS==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==8 .AND. LBEPKE/=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>0) CALL GROUND
      ENDIF
      if(flag.or.dbgloc) call banner(2,'SLBVST',131112 )
      NAMSUB= 'SLBVST'
      END
c-------------------------
      SUBROUTINE INIVST
      INCLUDE 'farray'
      INCLUDE 'qdata'
      INCLUDE 'satear'
      INCLUDE 'satgrd'
      INCLUDE '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',131112 )
C.... Turbulence-model constants:
      IF(.NOT.RSTM) THEN
        IF(CMU<=0.0)  CMU=0.5478
        IF(CD<=0.0)   CD= 0.1643
        IF(C1E<=0.0)  C1E=1.44
        IF(C2E<=0.0)  C2E=1.92
        IF(AK<=0.0)   AK= 0.41
        IF(EWAL<=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==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==2 .OR. IENUTA==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==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==10 .OR. IENUTA==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<0.0) THEN
        CALL ADDL1C(L0UTAU,IENUTA==5,'l0utau')
        IF(IENUTA==3 .OR. IENUTA==4 .OR. IENUTA==8 .OR.
     1     IENUTA==11) THEN
          NXNY=NX*NY
          IF(LBFONE==0) CALL GXMAKE(L0FONE,NXNY,'FONE')
          IF(LBFTWO==0) CALL GXMAKE(L0FTWO,NXNY,'FTWO')
          IF(LBFMU==0)  CALL GXMAKE(L0FMU,NXNY,'FMU ')
          IF(LBREYN==0.AND.IENUTA/=11) CALL GXMAKE(L0REYN,NXNY,'RE')
          IF(LBREYT==0.AND.IENUTA/=8) CALL GXMAKE(L0REYT,NXNY,'RET')
        ENDIF
      ENDIF
      if(flag) call banner(2,'INIVST',131112 )
      END
************************************************************************
      REAL FUNCTION ENUTPL(REYN)
      COMMON/TURB/CMU,CD,CMUCD,C1E,C2E,AK,EWAL
      COMMON/PLUSES/UPLUS,YPLUS  ! this common appears also in sknfg2
c
      RECEW=1.0/EWAL
      REYNHI=REYN*1.001         ! To stop search when within 0.1% of
      REYNLO=REYN*0.999         ! of specified Reynolds number
      IF(REYN<=1.0) THEN
        ENUTPL=0.0
        YPLUS=SQRT(REYN)
        UPLUS=YPLUS
        RETURN
      ELSEIF(REYN<=1.E1) THEN   ! this block of conditions is designed to
        UPSMALL=1.0             ! reduce the number of iterations.
        UPBIG=3.2               ! The values set for upbig and upsmall
      ELSEIF(REYN<=1.E2) THEN   ! are valid only for ak=0.41. ewal=8.6
        UPSMALL=3.1
        UPBIG=10.0
      ELSEIF(REYN<=1.E3) THEN
        UPSMALL=9.0
        UPBIG=16.0
      ELSEIF(REYN<=1.E4) THEN
        UPSMALL=15.0
        UPBIG=21.0
      ELSEIF(REYN<=1.E5) THEN
        UPSMALL=20.0
        UPBIG=26.0
      ELSEIF(REYN<=1.E6) THEN
        UPSMALL=25.
        UPBIG=31.
      ELSEIF(REYN<=1.E7) THEN
        UPSMALL=30.0
        UPBIG=36.0
      ELSEIF(REYN<=1.E8) THEN
        UPSMALL=35.0
        UPBIG=42.0
      ELSE
        UPSMALL=41.0
        UPBIG=100
      ENDIF
      DO I=1,100
        NUMIT=I
        UP=0.5*(UPBIG+UPSMALL)
        AKUP=AK*UP
        EXPKUP=EXP(AKUP)
        YP=UP+RECEW*(EXPKUP -( 1 + AKUP*(1 + AKUP*(0.5 +AKUP*(0.1666667
     1                         + 0.0444444*AKUP)))))
        RE=UP*YP
       IF(REREYNLO) THEN  ! but if the difference is small, don't bother
            GOTO 10
         ENDIF
         UPSMALL=UP
       ELSE                   !  up will need to decrease
         IF(RE