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#### MRM 19.01.17 Realisable k-e model 
C      COMMON/REALKE1/LBCMU,LBFIL1(10)
      COMMON/REALKE2/L0CMU,LBFIL2(10)   
C#### MRM 11.07.17 Menter k-w SST model         
      COMMON/KWMOD2/L0KWF1(7),L0BF1,L0BF2,L0BF3,L0KWF2(7)   
C   
C  List of turbulence models (according to IENUTA)
C  -----------------------------------------------       
c< laminar flow
c = 0   !< standard k-e model
c = 1   !< RNG high-Re k-e model 
c = 2   !< Chen-Kim high-Re k-e model 
c = 3   !< Lam-Bremhorst low-Re k-e model 
c = 4   !< Chen-Kim low-Re k-e model 
c = 6   !<  Simple eddy-viscosity multiplier for low Reynolds numbers 
c = 7   !< Two-scale high-Re k-e model 
c = 8   !< Two-layer low-Re k-e model 
c = 9   !< Another simple eddy-viscosity multiplier for low Reynolds numbers 
c = 10  !< Wilcox (1988) high-Re k-omega model 
c = 11  !< Wilcox (1988) low-Re k-omega model 
c = 12  !< MMK high-Re k-e model 
c = 13  !< Kato-Launder high-Re k-e model
c = 14  !< Realisable high-Re k-e model 
c = 15  !< Wilcox (2008) high-Re k-omega model 
c = 16  !< Wilcox (2008) low-Re k-omega model    
c = 17  !< Menter (1992) high-Re k-omega model 
c = 18  !< Menter (1992) low-Re k-omega model 
C = 19  !< high-Re k-omega SST model
C = 20  !< low-Re k-omega SST model 
c other turbulence models where IENUTA is not used:
c       !< LVEL model
c       !< model of constant effective viscosity      
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
C#### MRM 11.07.17 Menter k-w SST model             
        GOMEGA = AMAX1(F(L0OMEG+I),TINY)
        GKE    = F(L0KE+I)
C#### MRM 25.08.17 Wilcox (2008) k-w model 
        IF(IENUTA==15.OR.IENUTA==16) THEN  ! Wilcox (2008) k-w model
          CLIM   = 7./8.  
          GGEN1   = SQRT(F(L0GEN+I))  ! S=Sqrt(2.*SijSij) = sqrt(GEN1) 
          GOMLIM  = AMAX1(GOMEGA,CLIM*GGEN1/0.3)  
          KINVIST = GKE/GOMLIM 
        ELSE IF(IENUTA==19.OR.IENUTA==20) THEN  ! k-w SST model
          GGEN1   = SQRT(F(L0GEN+I))  ! S=Sqrt(2.*SijSij) = sqrt(GEN1) 
          GBF2    = F(L0BF2+I)
C          GBF3 = F(L0BF3+I)
C          GSF2 = GGEN1*GBF2*GBF3
          GSF2    = GGEN1*GBF2
          GOMLIM  = AMAX1(0.3*GOMEGA,GSF2)
          KINVIST = 0.3*GKE/GOMLIM
        ELSE                 ! other k-w models
          KINVIST= GKE/GOMEGA
        ENDIF
      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 '/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/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
C#### MRM 19.01.17 Realisable k-e model  
      COMMON/REALKE1/LBCMU,LBFIL1(10)
      COMMON/REALKE2/L0CMU,LBFIL2(10) 
C#### MRM 11.07.17 Menter k-w SST model 
      COMMON/KWMOD1/LBKWF1(7),LBBF1,LBBF2,LBBF3,LBSIGK,LBSIGW,LBKWF2(5)
      COMMON/KWMOD2/LBKWF3(7),L0BF1,L0BF2,L0BF3,L0SIGK,L0SIGW,LBKWF4(5)    
      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,
c#### dbs 31.10.10 kzxcy inserted. correctly?
     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',190117 )
        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(STORE(3)) L0UVW(3)=L0F(3)
        IF(STORE(5)) L0UVW(1)=L0F(5)
        IF(STORE(7)) L0UVW(5)=L0F(7)
        IF(.NOT.ONEPHS) THEN
          IF(STORE(4)) L0UVW2(3)=L0F(4)
          IF(STORE(6)) L0UVW2(1)=L0F(6)
          IF(STORE(8)) L0UVW2(5)=L0F(8)
        ENDIF  
      ENDIF  
      L0VISL  = L0F(VISL)
      IF(IGRND==3) THEN
        IF(STORE(12)) L0KE=L0F(12)
C#### JCL 20.10.16 need L0R2 if ENUTC>0 and STORE(R2)
        IF(ENUTC>0.0.AND.STORE(10)) L0R2=L0F(10)
C#### JCL 28.03.11 add missing settings for L0R1 and L0R2        
      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))
C#### JCL 03.08.12 get L0REYN for ENUT=GRND8        
        IF(LBREYN/=0) L0REYN= L0F(LBREYN)
C#### JCL 13.11.12 must define L0RATE and L0MIXL before call to KINVIST,
C####              definition in MFMSOLVER is far too late!        
      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)
C#### MRM 19.01.17 Realisable k-e model          
      ELSEIF(IENUTA==14) THEN
        IF(LBCMU/=0)  L0CMU = L0F(LBCMU)
C#### MRM 11.07.17 Menter k-w SST model    
      ELSEIF(IENUTA==19.OR.IENUTA==20) THEN
        IF(LBBF2/=0)  L0BF2 = L0F(LBBF2)   
C        IF(LBBF3/=0)  L0BF3 = L0F(LBBF3) 
      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
            GVIST = KINVIST(I)  
            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',190117 )
      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,YPLUSC,REWC 
     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)
C#### MRM 11.07.17 Menter k-w SST model       
      COMMON/KWMOD3/CWWALL,SIGK1,SIGK2,SIGW1,SIGW2,CWA1,CWA2,CWB1,CWB2,
C#### MRM 20.10.17 Revised k-w model - Cross-diffusion coefficient       
     1              BETAST,CDSIG
C#### MRM 03.02.17 Extend scalable wall function to k-f model 
C#### MRM 13.07.17 Menter k-w & k-w SST model        
      LOGICAL GETREALS,KEMOD,KFMOD,SOLVE,KFHIRE
      DIMENSION NAMPRP(30)
      EQUIVALENCE (NAMPRP(1),DENST1)
C.... Preliminaries:
      if(flag) call banner(1,'INIVST',120717 )
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#### JCL 10.01.17 get critical Y+ for scalable wall function
        IF(YPLUSC<=0.0) YPLUSC=11.126
C#### MRM 20.09.17 Turbulence-model coefficients clean-up for k-w models            
        KEMOD=SOLVE(12).AND.SOLVE(13)        
        KFMOD=SOLVE(12).AND.SOLVE(LBOMEG)
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)
          CALL GETSPD('KECONST','YPLUSC',1,YPLUSC,IDUM,.FALSE.,' ',IERR)
        ENDIF
        CMUCD=CMU*CD
C#### JCL 10.01.17 check validity of scalable wall function
C#### MRM/JCL 09.05.17 allow SCALWF for constant ENUT        
        IF(SCALWF.AND.ENUT/=GRND9) THEN
C#### MRM 20.09.17 Turbulence-model coefficients clean-up for k-w models                
c####          KEMOD=SOLVE(12).AND.SOLVE(13)
          IF(.NOT.KEMOD) THEN ! only for KE models
            SCALWF=.FALSE.
C#### MRM 03.02.17 Extend scalable wall function to k-f model 
C#### MRM 20.09.17 Turbulence-model coefficients clean-up for k-w models                
c####            KFMOD=SOLVE(12).AND.SOLVE(LBOMEG)
C#### MRM 13.07.17 Menter k-w & k-w SST model  
C#### MRM 25.08.17 Wilcox (2008) k-w model                
            KFHIRE=IENUTA==10.OR.IENUTA==15.OR.IENUTA==17.OR.IENUTA==19
            IF(KFMOD.AND.KFHIRE) SCALWF=.TRUE.! only high-RE
          ELSE 
            IF((IENUTA>2.AND.IENUTA<7).OR.(IENUTA>7.AND.IENUTA<10).OR.
     1                                                  IENUTA==11) THEN ! only high-RE
              SCALWF=.FALSE.
            ENDIF
          ENDIF
          IF(.NOT.SCALWF) THEN
            CALL WRIT40('Scalable wall function switched off because')
            CALL WRIT40('not compatible with current turbulence model')
          ENDIF
        ENDIF
        REWC=YPLUSC**2. ! critical wall Reynolds number
c#### dbs 11.03.10 compute quantities used in egwf
        RSCMCD=CMUCD**(-0.5) ; CDDAK2=2.0*CD/AK ; TAUDKE=SQRT(CMUCD)
        RTTDKE=SQRT(TAUDKE) ; AKC=AK*RTTDKE ; EWC=EWAL*RTTDKE        
C        
C#### MRM 30.10.17 Print out default model coefficients for IENUTA=3 & 8              
C.... Standard K-E model:
        IF(IENUTA==0) THEN
          IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                  ' Standard k-e model constants           ')
C#### MRM 30.10.17 Print out default model coefficients for IENUTA=3 & 8                
C.... Lam-Bremhorst low-Re K-E model:
        ELSEIF(IENUTA==3) THEN
            IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                  ' Lam-Bremhorst low-Re model constants   ')  
C.... Two-layer low-Re K-E model:
        ELSEIF(IENUTA==8) THEN
            IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                  ' Two-layer low-Re k-e model constants   ')  
        ELSEIF(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)
          IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                    ' RNG k-e turbulence model constants     ')
C#### MRM 30.10.17 Print out default model coefficients for IENUTA=3 & 8          
        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)
C#### MRM 30.10.17 Print out default model coefficients for IENUTA=3 & 8                
          IF(.NOT.NULLPR) THEN
            IF(IENUTA==2) CALL PR_TURCON(IENUTA,
     1                    ' Chen-Kim k-e turbulence model constants')
            IF(IENUTA==4) CALL PR_TURCON(IENUTA,
     1                    ' Chen-Kim low-Re k-e model constants    ')
          ENDIF  
        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)
          IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                    'Two-scale k-e model constants          ')       
        ELSEIF(IENUTA==10 .OR. IENUTA==11) THEN
C.... Wilcox (1988) high & low-Re k-w models
          C1E=5./9. ; C2E=0.075
          CWWALL=6.0 ! low-Re wall constant for w
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
     1                 '   ',IERR)
          KELIN= 0
C#### MRM 30.10.17 Improve printout of k-w model coefficients          
          IF(.NOT.NULLPR) THEN
            IF(IENUTA==10) CALL PR_TURCON(IENUTA,
     1                    ' Wilcox (1988) k-w model constants      ')
            IF(IENUTA==11) CALL PR_TURCON(IENUTA,
     1                    'Wilcox (1988) low-Re k-w model constants')            
          ENDIF  
        ELSEIF(IENUTA==14) THEN
C.... Realisable high-Re K-E model:
          C2E=1.9 
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
C#### MRM 15.02.17 Realisable model - force KELIN=3 for safety
C... KELIN=3 must be used because other KELIN options use a constant
C    CMU in the linearisation of the KE source/sink terms          
          KELIN=3   
          IF(.NOT.NULLPR) CALL PR_TURCON(IENUTA,
     1                    ' Realisable k-e model constants         ')
C#### MRM 25.08.17 Wilcox (2008) k-w model 
      ELSEIF(IENUTA==15 .OR. IENUTA==16) THEN
C.... Wilcox (2008) high & low-Re revised k-w models
          CDSIG=1./8. ! cross-diffusion source-term coefficient
          C1E=13./25. ; C2E=0.0708
          CWWALL=6.0 ! low-Re wall constant for w
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CDSIG',1,CDSIG,IDUM,.FALSE.,'   ',IERR)          
          CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
     1                 '   ',IERR)
          KELIN= 0   
C#### MRM 30.10.17 Improve printout of k-w model coefficients
          IF(.NOT.NULLPR) THEN
            IF(IENUTA==15) CALL PR_TURCON(IENUTA,
     1                    ' Wilcox (2008) k-w model constants      ')
            IF(IENUTA==16) CALL PR_TURCON(IENUTA,
     1                    'Wilcox (2008) low-Re k-w model constants')            
          ENDIF
C#### MRM 11.07.17 Menter (1992) Baseline k-w model               
      ELSEIF(IENUTA==17.OR.IENUTA==18) THEN
C.... Menter (1992) high & low-Re Baseline k-w models
          SIGK1 = 2.0 ; SIGW1 = 2.0   ! wall region limit
          SIGK2 = 1.0 ; SIGW2 = 1./0.856 ! outer region limit
          C1E=5./9. ; C2E=3./40.
          CWWALL=6.0 ! low-Re wall constant for w
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
     1                 '   ',IERR)
          BETAST=CMUCD
          CWB1=C2E ! = 0.075
          CWB2 = 0.0828 ; AK2 = AK*AK
          CWA1 = CWB1/BETAST - AK2/(SIGW1*SQRT(BETAST)) ! CWA1 = 0.5532=C1E
          CWA2 = CWB2/BETAST - AK2/(SIGW2*SQRT(BETAST)) ! CWA2 = 0.44
          KELIN= 0  
C#### MRM 30.10.17 Improve printout of k-w model coefficients          
          IF(.NOT.NULLPR) THEN
            IF(IENUTA==17) CALL PR_TURCON(IENUTA,
     1                    ' Menter Baseline k-w model constants    ')
            IF(IENUTA==18) CALL PR_TURCON(IENUTA,
     1                    'Menter Baseline low-Re k-w model        ')            
          ENDIF
C#### MRM 11.07.17 Menter k-w SST model               
        ELSEIF(IENUTA==19.OR.IENUTA==20) THEN
C.... Menter (1992) high & low-Re k-w SST models
C   Original Menter paper sets SIGK1=0.85 & SIGW1 = 0.65 ! check       
          SIGK1 = 2.0 ; SIGW1 = 2.0
          SIGK2 = 1.0 ; SIGW2 = 1./0.856 ! outer region limit          
          C1E=5./9. ; C2E=3./40.
          CWWALL=6.0 ! low-Re wall constant for w
          CALL GETSPD('KECONST','C1E',1,C1E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','C2E',1,C2E,IDUM,.FALSE.,'   ',IERR)
          CALL GETSPD('KECONST','CWWALL',1,CWWALL,IDUM,.FALSE.,
     1                 '   ',IERR)
          BETAST=CMUCD  ! CMUCD =0.09
          CWB1=C2E 
          CWB2 = 0.0828 ; AK2 = AK*AK
C### MRM 0.08.08.17 Error in CWA2 replace BETAW1 with CWB2          
          CWA1 = CWB1/BETAST - AK2/(SIGW1*SQRT(BETAST)) ! CWA1 = 0.5532
          CWA2 = CWB2/BETAST - AK2/(SIGW2*SQRT(BETAST)) ! CWA2 = 0.44
          KELIN= 0 
C#### MRM 30.10.17 Improve printout of k-w model coefficients          
          IF(.NOT.NULLPR) THEN
            IF(IENUTA==19) CALL PR_TURCON(IENUTA,
     1                    ' k-w SST turbulence model constants     ')
            IF(IENUTA==20) CALL PR_TURCON(IENUTA,
     1                    ' k-w SST low-Re model constants         ')
          ENDIF          
        ENDIF
C        IF(.NOT.NULLPR) THEN
C#### MRM 20/09.17 Turbulence-model coefficients clean-up for k-w models    
      ELSE
C#### JCL 10.01.17 check validity of scalable wall function
        IF(SCALWF) THEN
          SCALWF=.FALSE.
          CALL WRIT40('Scalable wall function switched off because')
          CALL WRIT40('not compatible with RSTM turbulence model')
        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',190117 )
      END
************************************************************************      
c#### dbs 19.06.12 New function for Spalding wall law
      REAL FUNCTION ENUTPL(REYN)                 
      COMMON/TURB/CMU,CD,CMUCD,C1E,C2E,AK,EWAL,YPLUSC,REWC
      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