c

C     SUBROUTINE GXWFUN provides "wall functions", of the kind
C     employed with turbulence models which do not solve the full
C     partial differential equations in a fine-grid region close to
C     a wall.
C
C.... This subroutine is called from GREX3 in the following locations:-
C     Group1 1, section 2; Group 13; and Group 19, sections 1 and 4.
C
C.... Five wall-function options are provided, namely:-
C
C     (a) Wall functions for fully turbulent flows according to
C         Blasius' Law (set COefficients to BLASIUS or GRND1);
C     (b) Log-law wall functions for smooth and rough walls (set
C         COefficients to GRND2, VALue to GRND2 for KE and EP);
C     (c) Generalised wall functions for smooth and rough walls
C         (set COefficients to GRND3, VALue to GRND3 for KE and EP).
C     (d) Fully-rough log-law wall functions (set COefficients to
C         GRND5, VALue to GRND5 for KE and EP);
C     (e) Spalding's wall function for smooth walls which is valid for
c         all Reynolds numbers
c
C     Options (b) and (c) allow the specification of a rough wall by
C     setting the PIL variable WALLA to a finite positive 'sand-grain'
C     roughness height. Option (d) defines a fully-rough wall in terms
C     of the effective roughness height. For information on how to
C     specify a different roughness height for each wall click on
C     'Roughness Wall Functions'in the PHOENICS Encyclopaedia.
C
C     A description of the wall-function formulae can be found in the
C     Encyclopaedia by clicking on 'Wall Functions'. Other
C     relevant Encyclopaedia entries are: WALL, WALLS, WALLCO, WALLA
C     and 'Roughness Wall Functions'.
C
      SUBROUTINE GXWFUN
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/satear'
      COMMON /GENI/  IGEN1(18),IPAT1,IGEN2(5),IPAT2,IGEN3(35)
      COMMON /UVWCOL/IUC1,IVC1,IWC1,IUF1(3),IUC2,IVC2,IWC2,IUF2(24)
      COMMON /SODAL/ DBG,LOGIC(8)
      COMMON/LRNTM1/ L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
     1               L0UD2,L0UD3,L0UD4
      COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
     1       /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
ccccc      COMMON /TURB3/ RTTDKE,AKC,EWC,ACON,TAUDKE 
     1       /LWFUN/CORGNR
C#### MRM 25.09.17 Introduce CWWALL for k-w models    
      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
      COMMON/TSKEMI/ JKP,JKT,JET,LBVOSQ,LBOMEG
     1      /TSKEMR/ GKTDKP
      COMMON /NAMFN/ NAMFUN,NAMSUB
C#### JCL/MRM 19.08.14 allow for 'd' constant in fully-rough log law
      LOGICAL DBG,LOGIC,CORGNR,LUVWC,DOSKIN,LBFCHV,MASFLOP,NEZ,GTZ,
     1        LPTDMN,EQZ,DOIT
      INTEGER COGRN,VALGRN
      CHARACTER*6 NAMFUN,NAMSUB,BUFF40*40
C
      NAMSUB='GXWFUN'
      NAMFUN='      '
      IF(IGR==1 .AND. ISC==2) THEN
C-----------------------------------------------------------------------
C---- Group 1 ---- Section 2 ---------------------------- Initialization
C
      ELSEIF(IGR==13) THEN
C-----------------------------------------------------------------------
C---- Group 13 ------------------------------------- Boundary conditions
        COGRN=ISC-1 ; VALGRN=ISC-12
        if(dbg) then
          call writ40('debug at entry to group 13 of gxwfun    ')
          call writ2i('indvar  ',indvar,'isc     ',isc)
          call writ2r('rttdke  ',rttdke,'taudke  ',taudke)
        endif
        IF(ITYPE==17.OR.ITYPE==18) THEN  ! E-W
          IDIR=3
        ELSEIF(ITYPE==19.OR.ITYPE==20) THEN ! N-S
          IDIR=1
        ELSEIF(ITYPE==21.OR.ITYPE==22) THEN ! H-L
          IDIR=5
        ENDIF
C.... Coefficients section
C     CO=BLASIUS (GRND1)  or LOGLAW (GRND3) or GENLAW (GRND3) or GRND5
c         or GRND6 for Spalding's wall law
        IF(COGRN>0.AND.COGRN<=6) THEN
C.... Compute skin-friction if not already done
          IF(DOSKIN(IPAT)) THEN
C.... Recalculate relative velocities for IPAT-patch, if
C     MB-FGE technique is in use.
            IF(CCM) THEN
              CALL UNRLVL(IPAT,IZSTEP)
            ELSEIF(GCV) THEN
              CALL BFRLVL(IPAT,IZSTEP)
            ENDIF
            CALL FNSKIN
          ENDIF
C.... Check for collocated velocity components:
          IF(CCM) THEN
            LUVWC=INDVAR==IUC1.OR.INDVAR==IVC1.OR.INDVAR==IWC1
     1          .OR.INDVAR==IUC2.OR.INDVAR==IVC2.OR.INDVAR==IWC2
          ELSEIF(GCV) THEN
            LUVWC=LBFCHV(INDVAR)
          ELSE
            LUVWC=.FALSE.
          ENDIF
          IF(INDVAR<9 .OR. LUVWC) THEN
C.... Compute coefficients for velocities
            CALL FNCOEF(1.0,1.0,.TRUE.)
C#### JCL 25.05.12 add LBOMEG and LBVOSQ
          ELSEIF(INDVAR>=14.AND.INDVAR/=JKP.AND.INDVAR/=JKT
     1           .AND.INDVAR/=JET.AND.INDVAR/=LBOMEG.AND.
     1                                       INDVAR/=LBVOSQ) THEN
C.... Compute coefficients for scalars
            CALL FNCOEF(PRNDTL(INDVAR),PRT(INDVAR),.FALSE.)
          ELSEIF(INDVAR==JET.OR.((INDVAR==JKP.OR.INDVAR==JKT).
     1           AND.COGRN==2)) THEN
            L0SU=L0F(LSU); L0AP=L0F(LAP)
            DO IX=IXF,IXL
              DO IY=IYF,IYL
                I=(IX-1)*NY+IY
                IF(.NOT.MASFLOP(I,IDIR).AND..NOT.MASFLOP(I,7)) THEN
                  F(L0SU+I)=0.0; F(L0AP+I)=1.1*FIXVAL
                ENDIF
              ENDDO
            ENDDO
          ELSEIF((INDVAR==12.OR.INDVAR==JKP.OR.INDVAR==JKT)
     1                                          .AND.COGRN==3) THEN
C.... Compute coefficients for ke or epsilon
            L0SU=L0F(LSU); L0CO=L0F(CO)
            DO IX=IXF,IXL
              DO IY=IYF,IYL
                I=(IX-1)*NY+IY
                IF(.NOT.MASFLOP(I,IDIR).AND..NOT.MASFLOP(I,7)) THEN
                  F(L0CO+I)=1.0
                ENDIF
              ENDDO
            ENDDO
          ENDIF
C.... Values section
        ELSEIF(VALGRN>=1) THEN
          IF(INDVAR==KE.OR.INDVAR==JKP.OR.INDVAR==JKT) THEN
            IF(VALGRN==2 .OR. VALGRN==5) THEN
C.... KE is here fixed to (wall-shear-stress)/(density*sqrt(cmucd))
              CALL FNPAXY(VAL,PVSTRS)
              IF(INDVAR==KE) THEN    
                CALL FN25(VAL,1.0/TAUDKE)
              ELSEIF(INDVAR==JKP) THEN
                CALL FN25(VAL,1.0/(TAUDKE*(1.+GKTDKP)))
              ELSEIF(INDVAR==JKT) THEN
                CALL FN25(VAL,GKTDKP/(TAUDKE*(1.+GKTDKP)))
              ENDIF
              L0VAR=L0F(INDVAR); L0VAL=L0F(VAL); L0CO=L0F(CO)
              RLX=-DTFALS(INDVAR)
              DO IX=IXF,IXL
                DO IY=IYF,IYL
                  I=(IX-1)*NY+IY
                  IF(MASFLOP(I,IDIR).OR.MASFLOP(I,7)) THEN
                    F(L0CO+I)=0.0; F(L0VAL+I)=0.0
                  ELSE
                    VAR=F(L0VAL+I)
                    IF(RLX>=0.0) VAR=(1.0-RLX)*F(L0VAR+I)+RLX*VAR
                    F(L0VAR+I)=AMIN1(VARMAX(INDVAR),AMAX1(VAR,
     1                               VARMIN(INDVAR))) ! equivalent to FN0 above
                  ENDIF
                ENDDO
              ENDDO
            ELSEIF(VALGRN==3) THEN
              CALL FNPAXY(VAL,PVGENR)
              CALL FN60(VAL,INDVAR)
              L0VAL=L0F(VAL); L0CO=L0F(CO)
              DO IX=IXF,IXL
                DO IY=IYF,IYL
                  I=(IX-1)*NY+IY
                  IF(MASFLOP(I,IDIR).OR.MASFLOP(I,7)) THEN
                    F(L0CO+I)=0.0; F(L0VAL+I)=0.0
                  ENDIF
                ENDDO
              ENDDO
            ENDIF
          ELSE
            IF(VALGRN==2.OR.VALGRN==3.OR.VALGRN==5) THEN
C.... Put reciprocal distance into EASP1
              CALL FNPAXY(EASP1,PVRCDS)
              IF(INDVAR==EP.OR.INDVAR==JET) THEN
C.... In this section, epsilon is set to cd*(ke)**1.5/(ak*waldis)
C     then compute value of EP for wall
C#### JCL/MRM 19.08.14 allow for 'd' constant in fully-rough log law
C.... In this section, epsilon is set to cd*(ke)**1.5/(ak*(waldis-d))
C####                IF(NEZ(WALLB)) THEN
                L0VAL=L0F(VAL); L0KE =L0F(KE); L0RECW =L0F(EASP1)
                L0VISL=L0F(VISL); L0STRS=L0PVAR(PVSTRS,IPAT,0); J=0
                DOIT=SCALWF.AND.INDVAR==EP.AND.EQZ(WALLA).AND.
     1                                                    VALGRN==2
C
                DO IX=IXF,IXL
                  DO IY=IYF,IYL
                    I=(IX-1)*NY+IY; J=J+1
                    DISTW=1./F(L0RECW+I)
C#### JCL/MRM 06.01.17 get wall distance consistent with scalable wall function
                    IF(DOIT) DISTW=MAX(DISTW, YPLUSC*F(L0VISL+I)/
     1                                           SQRT(F(L0STRS+J)+TINY))
C####                    IF(DOIT) then
C####                      distw0=distw
C####                      DISTW=MAX(DISTW, YPLUSC*F(L0VISL+I)/
C####     1                                           SQRT(F(L0STRS+J)+TINY))
C####                      if(distw.ne.distw0) then
C####                        write(14,'(''gxwfun '',4i4,1p2e12.3)') isweep,
C####     1                                         indvar,ix,iy,distw0,distw
C####                      endif
C####                    endif
                    IF(GTZ(WALLB)) THEN
                      DISWMD=AMAX1((DISTW-WALLB),2.0*WALLA)
                    ELSE
                      DISWMD=(DISTW-WALLB)
                    ENDIF
                    F(L0VAL+I)=CD*F(L0KE+I)**1.5/(AK*DISWMD+TINY) 
                  ENDDO
                ENDDO
C####                ELSE
C####                  CALL FN31(VAL,KE,EASP1,CD/AK,1.5,1.0)
C####                ENDIF                
              ELSEIF(INDVAR==LBVOSQ) THEN
C.... In this section, VOSQ is set to cd**2*ke/(ak*waldis))**2
C     then compute value of VOSQ for wall
C#### JCL 25.04.12 the comment says KE, but the arguments to FN31 say KE**2
C####              Change arguments so that KE is not squared.
C####                CALL FN31(VAL,KE,EASP1,(CD/AK)**2,2.0,2.0)
                CALL FN31(VAL,KE,EASP1,(CD/AK)**2,1.0,2.0)
              ELSEIF(INDVAR==LBOMEG) THEN
C#### MRM 13.07.17 Menter k-w SST model   
C#### MRM 05.07.17 Menter (1992) k-w model
C#### MRM 25.08.17 Wilcox (2008) k-w model
                IF(IENUTA==10.OR.IENUTA==15.OR.IENUTA==17
     1                       .OR.IENUTA==19) THEN
C.... In this section, OMEG=ke**0.5/(rttdke*ak*waldis)
C     where RTTDKE= CMUCD ** 0.25 
C####                  CALL FN31(VAL,KE,EASP1,1./(AK*RTTDKE),0.5,1.0)
C#### MRM 03.02.17 Extend scalable wall function to k-f model                     
                  L0VAL=L0F(VAL); L0KE =L0F(KE); L0RECW =L0F(EASP1)
                  L0VISL=L0F(VISL); L0STRS=L0PVAR(PVSTRS,IPAT,0); J=0
                  DOIT=SCALWF.AND.EQZ(WALLA)                                                            
                  DO IX=IXF,IXL
                    DO IY=IYF,IYL
                      I=(IX-1)*NY+IY; J=J+1
                      DISTW=1./F(L0RECW+I)
                      IF(DOIT) DISTW=MAX(DISTW, YPLUSC*F(L0VISL+I)/
     1                                           SQRT(F(L0STRS+J)+TINY))
C#### MRM 06.02.17 Allow for WALLB in k-f model
                      IF(GTZ(WALLB)) THEN
                        DISWMD=AMAX1((DISTW-WALLB),2.0*WALLA)
                      ELSE
                        DISWMD=(DISTW-WALLB)
                      ENDIF
                      F(L0VAL+I)=SQRT(F(L0KE+I))/(AK*RTTDKE*DISWMD+TINY)
                    ENDDO
                  ENDDO
C#### MRM 13.07.17 Menter k-w SST model
C#### MRM 25.08.17 Wilcox (2008) k-w model
                ELSEIF(IENUTA==11.OR.IENUTA==16.OR.IENUTA==18
     1                           .OR.IENUTA==20) THEN ! low-Re K-w models
C.... In this section, OMEG=CWWALL*visl/(c2e*waldis**2) with CWWALL=6 for 
C                      Wilcox k-w models & =10 for Menter k-w models
C#### MRM 25.09.17 Introduce CWWALL for k-w models  
                  CALL FN31(VAL,VISL,EASP1,CWWALL/C2E,1.0,2.0)
                ENDIF
              ENDIF
              L0VAR=L0F(INDVAR); L0VAL=L0F(VAL); L0CO=L0F(CO)
              RLX=-DTFALS(INDVAR)
              DO IX=IXF,IXL
                DO IY=IYF,IYL
                  I=(IX-1)*NY+IY
                  IF(MASFLOP(I,IDIR).OR.MASFLOP(I,7)) THEN
                    F(L0CO+I)=0.0; F(L0VAL+I)=0.0
                  ELSE
                    VAR=F(L0VAL+I)
                    IF(RLX>=0.0) VAR=(1.0-RLX)*F(L0VAR+I)+RLX*F(L0VAL+I)
                    F(L0VAR+I)=AMIN1(VARMAX(INDVAR),AMAX1(VAR,
     1                               VARMIN(INDVAR))) ! equivalent to FN0 above
                  ENDIF
                ENDDO
              ENDDO
            ELSEIF(INDVAR==14.OR.INDVAR==15) THEN
              IF(VALGRN==4) THEN
C.... Surface temperature set in TMP1/2 store for H1/2
                LBTMP1=LBNAME('TMP1'); LBTMP2=LBNAME('TMP2')
                LBTMP=LBNAME('T3')
                IF(LBTMP==0) LBTMP=ITWO(LBTMP1,LBTMP2,INDVAR==14)
                IF(LBTMP==0) THEN
                  IF(INDVAR==14) THEN
                    BUFF40='VAL=GRND4 for H1 requires STORE(TMP1)'
                  ELSE
                    BUFF40='VAL=GRND4 for H2 requires STORE(TMP2)'
                  ENDIF
                  CALL WRIT40(BUFF40)
                  CALL SET_ERR(240,'Error. See result file.',2)
                ENDIF
C... extract value from COVAL
                KPHCV=I0PHCV(LBTMP)-4
  60            KPHCV=KPHCV+4
                IF(KPHCV<=I0PHCV(LBTMP)+NPHCV(LBTMP)-4) THEN
                  IF(NINT(F(KPHCV+1))/=IPAT) GO TO 60
                  TWAL=F(KPHCV+3)+TEMP0 ! add TEMP0 as GASENT expects absolute value
                ENDIF
C... set enthalpy value
                KVPHI=L0F(VAL)
                DO IX=IXF,IXL
                  IADX=(IX-1)*NY
                  DO IY=IYF,IYL
                    I=IY+IADX
                    F(KVPHI+I)=GASENT(0,I,TWAL)
                  ENDDO
                ENDDO
                if(dbg) call prn('val ',val)
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ELSEIF(IGR==19 .AND. ISC==4) THEN
C-----------------------------------------------------------------------
C---- Group 19 ---- Section 4 ----------------------- Start of iteration
        DO IWL=1,NMWALL
          IPAT=NINT(ABS(F(L0WALL+IWL)))
          IF(LPTDMN(IPAT)) CALL SETSKN(IPAT)
        ENDDO
        IF(CORGNR .AND. (ISWEEP>FSWEEP.OR.ITHYD>1)) THEN
C.... The in-slab store for the strain-rate squared (set in subroutine
C     GXGENK) is overwritten here with what is set in Group 13 above
C     (subroutine FNSKIN), so as to ensure that the wall generation
C     function is correct in the near-wall reqions.
C.... CORGNR is set in EARTH and it is TRUE, if GRND1,2,3 wall-functions
C     are activated for at least one wall-patch.
C.... The correction procedure relies on the fact that PVDISS storage is
C     zeroised at the beginning of the run and should not be modified
C     except by FNSKIN.
          L0GEN1=L0F(GEN1)
          DO IWL=1,NMWALL
            IPAT=NINT(ABS(F(L0WALL+IWL)))
            IF(.NOT.LPTDMN(IPAT)) CYCLE
            L0G=L0PVAR(PVDISS,IPAT,0)
            CALL GETPAT(IPAT,IDUM,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2)
            IF(L0G/=0 .AND. IZSTEP>=IZ1 .AND. IZSTEP<=IZ2) THEN
              IPLUS=(IX1-2)*NY ; J=0
              DO IX=IX1,IX2
                IPLUS=IPLUS+NY
                DO IY=IY1,IY2
                  I=IY+IPLUS; J=J+1
                  F(L0GEN1+I)=F(L0G+J)
                ENDDO  
              ENDDO  
            ENDIF
          ENDDO  
        ENDIF
      ENDIF
      NAMSUB='gxwfun'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... FNSKIN calculates the dimensionless skin-friction factor
C     SKINF=stress/(rho * relvel**2) from one of:-
C     the Blasius law,         with COVAL arguments BLASIUS or GRND1;
C     the logarithmic law,     with COVAL arguments BLASIUS or GRND1;
C     the generalised log law, with COVAL arguments BLASIUS or GRND1.
c
      SUBROUTINE FNSKIN
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/lbnamer'
      COMMON /SODAL/DBG,LOGIC(8)
      COMMON /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
     1       /LWFUN/CORGNR
      COMMON /INDPV/NPVMX,NIMAX,L0PVS,L0PV(30)
      COMMON /GENL/ LGF1(56),TURB,LGF2(3)
      COMMON /GENI/LFG1(34),ISKIN,ISTAN,IYPLS,ISTRS,LFG2(22)
      COMMON /NAMFN/NAMFUN,NAMSUB
      LOGICAL CORGNR,GTZ,ROUGH,LGF1,TURB,LGF2,DBG,LOGIC,sld
      INTEGER COGRN
      CHARACTER*6 NAMFUN,NAMSUB
      COMMON /IEGWF/LBHTCO,LBHFLX,LBSHRX,LBSHRY,LBSHRZ,LBDIRX,LBDIRY,
     1             LBDIRZ
      LOGICAL EQZ
C
      NAMSUB='FNSKIN'
      COGRN=ISC-1
      if(dbg) then
        call writ40('Entry to FNSKIN                         ')
        call writ2i('INDVAR  ',indvar,'COGRN   ',cogrn)
        call prnpat('Rel Vel ',pvrlvl,ipat)
      endif
C.... LD11 contains laminar viscosity divided by wall distance,
C     i.e. the laminar wall coefficient.
      L0MUDD=L0F(LD11) ; L0SK=L0PVAR(PVSKIN,IPAT,0)
      L0RLVL=L0PVAR(PVRLVL,IPAT,0) ; L0STRS=L0PVAR(PVSTRS,IPAT,0)
      L0RCDS=L0PVAR(PVRCDS,IPAT,0); L0GENW=0
      IF(CORGNR) THEN
        L0GEN1=L0F(GEN1) ; L0GENR=L0PVAR(PVGENR,IPAT,0)
        L0GENW=L0PVAR(PVDISS,IPAT,0) ; IF(ENUT<0.0) L0VIST=L0F(VIST)
      ENDIF
C.... Logarithmic law options (GRND2 or GRND3):
      L0YPLS=0 ; L0ROGH=0
      IF(L0PV(PVYPLS)/=0) L0YPLS=L0PVAR(PVYPLS,IPAT,0)
      IF(L0PV(PVROGH)/=0) L0ROGH=L0PVAR(PVROGH,IPAT,0)
C#### JCL 23.04.10 bring in ROUG from PHI file as in Windsim ground, but
C####              allow patchwise store to overwrite
      IF(ISROUG.AND.L0ROGH==0) L0ROGH=-L0F(LBROUG)
      ROUGH=GTZ(WALLA) .OR. L0ROGH/=0
      SANDRG=WALLA
      if(dbg) call writ1l('ROUGH   ',rough)
C... switch to GRND2 wall function if no KE
      IF(COGRN==3.AND..NOT.STORE(KE)) COGRN=2
      IF(COGRN==2) THEN
C....  is possible to test the effect of the number of iterations
C      the solution by setting ISKINA in SATELLITE
        NUMITS=5+ISKINA
C####        SHALFM=1.0/11.5 ! never used
      ELSEIF(COGRN==3) THEN
        L0KE=L0F(KE)
      ENDIF
C.... Start of DO loop
      J=0; IPLUS=(IXF-2)*NY
      DO 20 IX=IXF,IXL
        IPLUS=IPLUS+NY
        DO 20 IY=IYF,IYL
          I=IY+IPLUS; J=J +1
C.... Reynolds number=relvel/(kinematic viscosity / wall distance)
          RELVEL=F(L0RLVL+J);  VISDDIS=F(L0MUDD+I)
          IF(VISDDIS<=0.0) THEN
            SKINF=0.0; GO TO 20
          ENDIF                      
          REYN=RELVEL/(VISDDIS + TINY)
          SKINFL=1.0/(REYN + 1.E-20)
C####          SLHALF=SQRT(SKINFL) ! never used
          IF(COGRN==1) THEN
            SKINF=BLASIUS(REYN)                ! Blasius law (GRND1)
          ELSEIF(COGRN==2) THEN
C.... Logarithmic law (GRND2, smooth or rough)
C#### JCL/MRM 06.01.17 apply scalable wall function (i.e. Y+>=11.126) to prevent
C####                  silly values from cells with very small distw
            IF(SCALWF.AND..NOT.ROUGH) REYN=MAX(REYN, REWC)
C####            IF(SCALWF.AND..NOT.ROUGH) then
C####              reyn0=reyn
C####              REYN=MAX(REYN, REWC)
C####              if(reyn/=reyn0) then
C####                write(14,'(''fnskin '',4i4,1p2e12.3)') isweep,indvar,
C####     1                                              ix,iy,reyn0,reyn
C####              endif 
C####            endif
            SKINFL=1.0/(REYN + 1.E-20)
            IF(ROUGH) THEN
              IF(L0ROGH/=0) THEN    ! roughness stored
                IF(L0ROGH>0) THEN  ! in patchwise store
                  SANDRG=F(L0ROGH+J)
                ELSE                  ! in 3D store via lbnamer
                  SANDRG=F(-L0ROGH+I)
                ENDIF
              ENDIF
              REYROU=REYN*SANDRG*F(L0RCDS+J)
            ELSE
              EREYNO=EWAL*REYN
            ENDIF
C... the function SKNFG2 returns max(Sturb, Slam)            
            SKINF = SKNFG2(REYN,NUMITS,ISKINB,ROUGH,REYROU,EEFF,SKINT)
            IF(L0YPLS/=0) F(L0YPLS+J)=REYN*SQRT(SKINF)
C... save E into EASP2 for use in Stanton Number            
            IF(ROUGH) THEN
              L0EA2=L0F(EASP2)
              F(L0EA2+I)=EW(SQRT(SKINT)*REYROU)
            ENDIF
            SKINF=SKINT ! return Sturb
          ELSEIF(COGRN==3) THEN
C.... Generalised logarithmic law (GRND3, smooth or rough)
            SQRK=SQRT(F(L0KE+I))
            SQKDVS=SQRK*F(L0RCDS+J)/F(L0MUDD+I)
            IF(ROUGH) THEN
              RERDH=RTTDKE*SQKDVS
              IF(L0ROGH/=0) THEN    ! roughness stored
                IF(L0ROGH>0) THEN  ! in patchwise store
                  SANDRG=F(L0ROGH+J)
                ELSE                  ! in 3D store via lbnamer
                  SANDRG=F(-L0ROGH+I)
                ENDIF
              ENDIF
              TERM=EW(RERDH*SANDRG)*RERDH
            ELSE
              TERM=EWC * SQKDVS
            ENDIF
            SKINF=AKC*SQRK/(ALOG(1.01+TERM/F(L0RCDS+J)) * F(L0RLVL+J))
            IF(L0YPLS/=0) F(L0YPLS+J)= REYN*SQRT(SKINF)
          ELSEIF(COGRN==5) THEN
C... Fully-rough logarithmic law for atmospheric b.l's (GRND5)
            IF(L0ROGH/=0) THEN    ! roughness stored
              IF(L0ROGH>0) THEN  ! in patchwise store
                SANDRG=F(L0ROGH+J)
              ELSE                  ! in 3D store via lbnamer
                SANDRG=F(-L0ROGH+I)
              ENDIF
            ENDIF
C#### MRM 29.11.11 Protect against wall distance <= roughness height
C... Zmin = Zo*exp{kappa/(i*(CmuCD)^0.5)} where i=1.0: - to prevent large skinf,ke and ep
C#### JCL 18.12.12 add tiny in case SANDRG=0
C#### JCL/MRM 19.08.14 allow for 'd' constant in fully-rough log law
            IF(EQZ(WALLB)) THEN
              ARG = AMAX1(1./(F(L0RCDS+J)*SANDRG+TINY),2.0)
            ELSE
              DISTW=1./F(L0RCDS+J)
              IF(LTZ(WALLB)) THEN
                ARG=(DISTW-WALLB)/(SANDRG+TINY)
              ELSE
                ARG=AMAX1((DISTW-WALLB)/(SANDRG+TINY),2.0)
              ENDIF
            ENDIF  
            SHALF=AK/(ALOG(ARG))
            SKINF=SHALF*SHALF
            IF(L0YPLS/=0) F(L0YPLS+J)= REYN*SQRT(SKINF)
          ELSEIF(COGRN==6)  THEN
C... Spalding's wall law
C... function SKNFG2 returns max(Sturb, Slam)            
            AAA=SKNFG2(REYN,5+ISKINA,ISKINB,.FALSE.,REYN,EEFF,SKINT)
            SKINF=SKINT
          ENDIF
          F(L0SK+J)=AMAX1(SKINF, SKINFL)
C.... Put into patchwise store PVSTRS the turbulent shear stress
C     divided by density
          IF(F(L0RLVL+J)>=1.E4) THEN ! If relative velocity > 1.e4 set to zero'
            F(L0RLVL+J)=0
          ENDIF  
          F(L0STRS+J)=F(L0SK+J)*F(L0RLVL+J)**2
          IF(CORGNR .AND. L0GENW/=0) THEN
C.... Here the generation rate computed in GREX is modified. The
C     generation rate to be used in the K-balance equation is placed in
C     the patch-wise store PVGENR after subtraction of the dissipation
C     rate which is calculated but not stored.
            GENRTN=F(L0STRS+J)*F(L0RLVL+J)
C.... The velocity-squared term is adjusted here to an appropriate
C     value for use in the enthalpy equation
            IF(SLD(I)) THEN ! Ensure that GENK is exactly zero in solid cells 
              VLGRD2=0.0    ! to avoid spurious TEM1/H1 sources .       
            ELSE
              VLGRD2=GENRTN*0.5*F(L0RCDS+J)
            ENDIF
            IF(ENUT>0.0) THEN
              F(L0GEN1+I)=VLGRD2/ENUT
            ELSE
C... In equilibrium it seems gen1 = Urel*Utau/(2.*kappa*disw^2)                
              F(L0GEN1+I)=VLGRD2/(F(L0VIST+I)+TINY)
            ENDIF
            F(L0GENW+J)=F(L0GEN1+I)
            IF(COGRN==3) THEN
C.... Calculate the dissipation rate from a formula which makes
C     DISS equal GENRTN then tau/rho equals the constant TAUDKE. Note that:
C        GENRTN=(tau/rho) * relvel
C        DISS=(tau/rho)**1.5/sqrt(skinf)
              DISS=(F(L0KE+I)*TAUDKE)**1.5*(ALOG(1.01+TERM/
     1               F(L0RCDS+J)))/AK
              F(L0GENR+J)=GENRTN - DISS
            ENDIF
          ENDIF
   20 CONTINUE
      IF(ISTRS/=0) CALL FNPAXY(ISTRS,PVSTRS)
      IF(ISKIN/=0) CALL FNPAXY(ISKIN,PVSKIN)
      IF(IYPLS/=0) CALL FNPAXY(IYPLS,PVYPLS)
!      ishr=0
!      if(lbshrx>0.and.indvar==3) then
!        ishr=lbshrx
!      elseif(lbshry>0.and.indvar==5) then
!        ishr=lbshry
!      elseif(lbshrz>0.and.indvar==7) then
!        ishr=lbshrz
!      endif
!      if(ishr>0) then
!        call fnpaxy(ishr,pvstrs)
!        call mupaxy(ishr,pvrlvl)
!        call fn55(ishr,den1,ishr,indvar,1.0)
!      endif  
C.... Set DOSKIN to false to prevent several calls to FNSKIN at the
C     same iteration   (sweep ?)
      CALL CLRSKN(IPAT)        ! does what?
      if(dbg) then
        call writ40('Exit from FNSKIN                        ')
        call prnpat('SkinFrCo',pvskin,ipat)
        call prnpat('TAU/RHO ',pvstrs,ipat)
        if(corgnr) then
          call prnpat('LGENWall',pvdiss,ipat)
          if(cogrn==3) call prnpat('GEN-DISS',pvgenr,ipat)
        endif
      endif
      NAMFUN='fnskin'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... FNCOEF supplies the wall coefficient/density, for smooth and rough
C     walls, and for both velocities and scalars.
c
      SUBROUTINE FNCOEF(PRLAM,PRTRB,VELVAR)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/lbnamer'
      INCLUDE  '/phoenics/d_includ/parear'
      COMMON /GENI/IGEN1(34),ISKIN,ISTAN,IYPLS,ISTRS,IGEN2(11),ITEM1,
     1             ITEM2,IGEN3(9)
      COMMON /GENL/LGN1(34),VELCTY,UVLCTY,VVLCTY,WVLCTY,LGN2(22)
      COMMON /SODAL/DBG,LSD1(5),STGEOM,LSD2(2)
      COMMON /RDATA/TINY,RSGFIL(84) /RDA3/PRNDTL(150)
      COMMON /IDATA/NX,NY,NZ,IDFIL1(66),I0PAT,NUMPAT,IDFIL2(49)
      COMMON /INDPV/NPVMX,NIMAX,L0PVS,L0PV(30)
      COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
     1       /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
cccc      COMMON /TURB3/RTTDKE,AKC,EWC,ACON,TAUDKE
      COMMON /NAMFN/NAMFUN,NAMSUB
      COMMON /RSTMCM/L0GTRS,L0PRTR
      COMMON /LDATA/ LDAT1(7),XCYCLE,LDAT2(23),NULLPR,LDAT3(52)
      LOGICAL LDAT1,LDAT2,LDAT3
      COMMON/IEGWF/LBHTCO,IEGFIL(7)
      LOGICAL LGN1,VELCTY,UVLCTY,VVLCTY,WVLCTY,LGN2,DBG,LSD1,STGEOM,
     1        LSD2,QEQ,PRSTAN,GTZ,ROUGH,VELVAR,SLD,GRN,GETPRL,LDOSTG,
     1        LSOLID,LSLDR,XCYCLE,PRHTCO,NULLPR
      INTEGER COGRN
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMFUN='FNCOEF'
      L0CO=L0F(CO)
C.... On entry, CO contains the laminar wall coefficient/density.
C     On exit, it contains the larger of this and the turbulent wall
C     coefficient/density, that is to say the Stanton number times
C     the relative velocity. Multiplication by density is effected in EARTH
c
      L0SK=L0PVAR(PVSKIN,IPAT,0)
      L0RLVL=L0PVAR(PVRLVL,IPAT,0)
      IF(VELVAR) THEN
C.... Coefficient for velocities
        if(dbg) then
          call writ40('Entry to FNCOEF from GXWFUN for velocity')
          call writ1i('INDVAR  ',indvar)
        endif
        LDOSTG=.FALSE.
        IF(STGEOM .AND. VELCTY) THEN
C.... Use averaging of the coefficient for staggered velocities,
C     when "there is no correspondent velocity at the patch edge":
          IJ=IYL+(IXL-1)*NY
          IF(UVLCTY) THEN
            LDOSTG=.NOT.XCYCLE.AND.(IXL==NX .OR.
     1               (IXL0) THEN  ! in patchwise store
                        SANDRG=F(L0ROGH+J)
                      ELSE                  ! in 3D store via lbnamer
                        SANDRG=F(-L0ROGH+I)
                      ENDIF
                    ENDIF
                    RTTDKH=RTTDKE*SANDRG
                    EWR=EW(RTTDKH*SQRK*F(L0RCDS+J)/F(L0MUDD+I))
                    FACTOR=CONST*PFUNR(EWR,PJAY)
                  ELSE
                    FACTOR=CONST*PJAY
                  ENDIF
                  RSTT=FACTOR*SKINF*F(L0RLVL+J)/SQRK
                ENDIF
C#### JCL/MRM 14.01.14 add TINY in case RSTT is exactly -1 leading to STAN=Inf then
C####                  CO=NaN e.g Grontmij bucket case US/12/021/13        
                STAN=SKINF/(PRTURB*(1.0 + RSTT)+TINY)
              ENDIF
              F(L0CO+I)=AMAX1(F(L0CO+I), STAN*F(L0RLVL+J))
              IF((L0STAN/=0.AND.PRSTAN).OR.PRHTCO) THEN
                STNLOC=AMAX1(STAN, F(L0CO+I)/(F(L0RLVL+J)+TINY))
                IF(L0STAN/=0.AND.PRSTAN) F(L0STAN+J)=STNLOC
                IF(PRHTCO) F(L0HTCO+I)=STNLOC*F(L0RLVL+J)
     1                                  *F(L0RHO1+I)*F(L0CP1+I)
              ENDIF
            ENDIF
          ENDDO
        ENDDO
        IF(GETPRL) PRNDTL(INDVAR)=PRLSTR
      ENDIF
      IF(ISTAN/=0) THEN
        CALL FNPAXY(ISTAN,PVSTAN)
        CALL FN22(ISTAN,1.E-10)
      ENDIF
      if(dbg) then
        call writ40('Exit from FNCOEF called from GXWFUN     ')
        call prn('COEF',co)
      endif
      NAMFUN='fncoef'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... This function calculates EW for rough walls of a certain class
      FUNCTION EW(RER)
      COMMON/TURB/TRBFIL(6),EWAL,YPLUSC,REWC
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      NAMFUN='EW    '
      EW=EWAL
      IF(RER>=3.7) THEN
        EW=29.7/RER
        IF(RER<=100.0) THEN
          X=0.02248*(100.0 - RER)*RER**(-0.564)
          ALF=1.0+X*X*(2.0*X-3.0)
          EW=EW/SQRT(ALF + (1.0 - ALF)*(EW/EWAL)**2)
        ENDIF
      ENDIF
      NAMFUN='ew    '
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... This function calculates the sublayer-resistance factor for
C     rough walls.
      FUNCTION PFUNR(EWR,PFUNS)
      COMMON/TURB/TRBFIL(6),EWAL,YPLUSC,REWC
      COMMON/RDA3/PRNDTL(150)
      COMMON/IGE/IXF,IXL,IYF,IYL,IGESP1(9),INDVAR,IGESP2(11)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      NAMFUN='PFUNR '
      EWDEWR=EWAL/EWR
      PRLPWR=PRNDTL(INDVAR)**0.695
      PFUNR=3.15*PRLPWR*(ABS(EWDEWR - 1.)/EWAL)**0.359 +
     1                                PFUNS*EWDEWR**(-0.6)
      NAMFUN='pfunr '
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXYPLS prints-out the near-wall values of Y+; the local shear
C     stress divided by the local density; the local skin friction
C     factors S; and, the local Stanton numbers for all surfaces defined
C     by wall-type patches, except those whose names begin with NF.
C.... GXYPLS is called from GREX3 in Group 19 Section 7, when YPLS (or
C     WALPRN) is set to T in the Q1 file.
c
      SUBROUTINE GXYPLS
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
      COMMON /GENI/IG1(18),IPAT1,IG2(5),IPAT2,IG3(9),ISKIN,ISTAN,IYPLS,
     1             ISTRS,IG4(11),ITEM1,ITEM2,IG5(9)
     1       /GENL/LGF1(15),SWPRIN,LGF2(44)
     1       /NAMFN/NAMFUN,NAMSUB
      LOGICAL LGF1,SWPRIN,LGF2
      CHARACTER NAMFUN*6,NAMSUB*6
C
      NAMSUB='GXYPLS'
C.... Print out Y+ etc. at last sweep
      IF(NULLPR) RETURN
      IF((ISWEEP==LSWEEP .OR. ENUFSW) .AND. SWPRIN) THEN
        CALL WRITBL
        CALL WRITST
C.... Don't print-out, if 3D-storage is provided
        IF(IYPLS==0 .OR. ISTRS==0 .OR. ISKIN==0 .OR. ISTAN==0)
     1    CALL WRIT40(' Patch-wise Printout from GXYPLS        ')
        IF(IYPLS==0) CALL PRNPAT('Y+      ',PVYPLS,0)
        IF(ISTRS==0) CALL PRNPAT('Tau/Rho ',PVSTRS,0)
        IF(ISKIN==0) CALL PRNPAT('Sk.Fr.Co',PVSKIN,0)
        IF(IVPRST/=0 .OR. SOLVE(H1) .OR. SOLVE(ITEM1) .OR.
     1                      SOLVE(H2) .OR. SOLVE(ITEM2)     ) THEN
          DO 10 IWAL=1,NMWALL
            IWPT=NINT(ABS(F(L0WALL+IWAL)))
            IF(NAMPAT(IWPT)/='NF' .AND. ISTAN==0)
     1        CALL PRNPAT('STANTON ',PVSTAN,IWPT)
  10      CONTINUE
        ENDIF
        CALL WRITST
        CALL WRITBL
      ENDIF
      NAMSUB='gxypls'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      REAL FUNCTION SKNFG2(REYNO,NIT,ISKINB,ROUGH,RER,EEFF,SKINT)
C.... This function calculates the skin-friction factor for the GRND2
C     option. It is normally called from within EARTH, when EGWF=T, but
C     it may also be called from GROUND
      INCLUDE '/phoenics/d_includ/farray'
      COMMON /TURB/ CMU,CD,CMUCD,C1E,C2E,AK,EWAL,YPLUSC,REWC
     1       /LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
c      COMMON TURB3/RTTDKE,AKC,EWC,ACON,TAUDKE
      COMMON/PLUSES/UPLUS,YPLUS  ! link with function enutpl
      LOGICAL ROUGH
C
C.... Preliminaries.
      IT=0 ; EEFF=EWAL ; SH=0.05 ; RAK=1./AK ; EREYN=EWAL*REYNO
      IF(ISKINB/=0) THEN  ! where set?
C.... The following loop solves the Spalding interpolating wall-law
C     equation for the Ku+, and hence the skin-friction factor
c#### dbs 19.06.12 There is much unexplained here, including:
c       1. The justification for the roughness treatment, for which no
c          no provison was made in the Spalding document.
c       2. The appearance of the REYNO>400. condition.
c       3. The disappaearance of the k_uplus**4 term which is present in 
c          the encyclopaedia article 
C       I have deactivated it, replacing its function by a call to ENUTPL
c       but leaving the coding in place in case back-tracking is needed
        GO TO 1000         
        AK2RE=AK*AK*REYNO ; EWDAK=EWAL*RAK ; AKDEW=1./EWDAK
        EWAKRE=AK*EREYN
        IF(ROUGH) AKRER=AK*RER
        AKUP=AK/SH
  10    IT=IT+1   ! start new iteration
        AKUPL=AKUP ; RAKUP=1./AKUP
        IF(ROUGH) THEN
          EEFF=EW(AKRER*RAKUP) ; EWDAK=EEFF*RAK ; AKDEW=1./EWDAK
          EWAKRE=EEFF*AK*REYNO
        ENDIF
        AKUP2=AKUP*AKUP
        IF(REYNO>400.0) THEN
          AKUP=(1. + LOG(1.01 + EWAKRE*RAKUP))/(1. + RAKUP)
        ELSE
          ARG=((0.1666667*AKUP + 0.5)*AKUP + 1.)*AKUP + 1.
          GX=(EXP(AKUP) - ARG)*AKDEW
          AKUP=((GX + 1.)*AKUP2 + AK2RE) / (2.*AKUP+(1. + AKUP)*GX +
     1           0.04166667*AKDEW*AKUP2**2)
        ENDIF
        IF(IT==NIT) THEN
          GO TO 20
        ELSEIF(IT==1 .OR. ABS(AKUP/AKUPL-1.)>1.E-4) THEN
          GO TO 10
        ENDIF
  20    SH=AK/AKUP
c#### dbs 19.06.12 SH appears to be sqrt(shear stress/density)/velocity
c                  i.e. the reciprocal of uplus
 1000 enutplus=enutpl(reyno)
c#### dbs 19.06.12 Tests show that the new coding and the old are in 
c                  satisfactory agreement
cccc        write(14,*) 'reyno,sh,uplus,sh*uplus=', reyno,sh,uplus,sh*uplus
        SH=1.0/UPLUS  
      ELSE
C.... The following loop solves the log-law equation for SQRT(SKINF)
  30    IT=IT+1
        SHL=SH
C#### JCL 23.04.12 pass EEFF back out so can be used to get Stanton Number        
C####        IF(ROUGH) EREYN=EW(SH*RER)*REYNO
        IF(ROUGH) THEN
          EEFF=EW(SH*RER); EREYN=EEFF*REYNO
        ENDIF  
        SH=AK/ALOG(1.01 + EREYN*SH)
        IF(IT==NIT) THEN
          GO TO 40
        ELSEIF(IT==1 .OR. ABS(SH/SHL-1.)>1.E-4) THEN
          GO TO 30
        ENDIF
      ENDIF
  40  SKINF=SH*SH
      SKINT=SKINF
      SKNFG2=AMAX1(SKINF,1./REYNO)
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
      FUNCTION STNFG2(REYNO,SKINF,PRL,PRT,PFAC,ROUGH,EEFF)
      INCLUDE '/phoenics/d_includ/grdear'
      COMMON /RDATA/TINY,RDFIL1(84)
C.... This function calculates the Stanton No. for the GRND2 option. 
      COMMON /RDA3/PRNDTL(150)
      LOGICAL ROUGH
C
      IF(ROUGH) THEN
        PRLS=PRNDTL(INDVAR)   ! save Q1 value of PRNDTL()
        PRNDTL(INDVAR)=PRL    ! set to current (local) value
        PEFF=PFUNR(EEFF,PFAC) ! PFUNR uses PRNDTL() in common as Prandlt Number
        PRNDTL(INDVAR)=PRLS   ! restore Q1 value
      ELSE                  
        PEFF=PFAC
      ENDIF
C#### JCL 16.01.14 add TINY to prevent zero-divide error when PEFF*SQRT(SKINF)
C####              is exactly -1. (e.g. US 01/015/14)      
      STAN=SKINF/(PRT*(1.0 + PEFF*SQRT(SKINF))+TINY)
C#### JCL 16.02.16 add TINY to prevent Infinity Stanton Number when REYNO and 
C####              PRL are both extremely small (e.g. US/01/020/16)      
      STNFG2=AMAX1(STAN,1./(REYNO*PRL+TINY))
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      REAL FUNCTION BLASIUS(REYN)    ! returns  coeff=shear stress / 
      COMMON /RDATA/TINY,RDFIL1(84)  !          (density  velocity ** 2)
      ! according to Blasius's smooth-wall law:  coeff=0.0395*REYN**(-0.25)
c!!!! dbs 09.03.10 There is no point in re-calculating the critical REYN
c                  on  each entry, here and in CHILTON
      CTURB=0.0395
cc      IF(REYN>CTURB**(-1.333)) THEN
      IF(REYN>74.366) THEN
        BLASIUS=CTURB * REYN ** (-0.25)
      ELSE
        BLASIUS=1.0/(REYN + TINY)  ! low-Reynolds number over-ride
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      REAL FUNCTION CHILTON(REYN,PRANDTL) ! returns the Stanton number 
      COMMON /RDATA/TINY,RDFIL1(84)       ! according to the Chilton-Colburn 
                                          ! heat/mass-transfer formula
      CTURB=0.0395
      IF(REYN>74.366) THEN
        BLASIUS=CTURB * REYN ** (-0.25)
        CHILTON=BLASIUS * PRANDTL ** (-0.6667)
      ELSE
        BLASIUS=1.0/(REYN + TINY)
        CHILTON=BLASIUS / PRANDTL  ! low-Reynolds number over-ride
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c