cgx2phs.htm c c

c  -- file name gx2phs.htm   030100
C     SUBROUTINE GXDISP is called from Group 13, only if the
C     PATCH name starts with KEDI.  The PATCH type should be CELL,
C     and COeff = GRND1, VAL = 0 for both KE and EP for the Mostafa
C     & Mongia model; COeff = GRND2, VAL=0.0 for the Chen and Wood
C     model; and COeff = FIXFLU, VAL = GRND3 for the bubble-induced
C     turbulence production model ( for more details see PHENC entry
C     'TURBULENCE MODELLING FOR TWO-PHASE FLOWS').
c
      SUBROUTINE GXDISP
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      common/genl/dbfil1(14),debgtz,dbfil2(45)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      EQUIVALENCE (JPRODB,EASP3)
      LOGICAL EQZ
      logical dbfil1,debgtz,dbfil2
C
      IF(IXL.LT.0) RETURN
      NAMSUB='GXDISP'
      IF(INDVAR.EQ.KE.OR.INDVAR.EQ.EP) THEN
        IF(ISC.EQ.2.OR.ISC.EQ.3) THEN
C.... Put large-scale turbulent motion timescale Te into EASP2
          CONST=0.35
          IF(IGR.EQ.3) CONST=0.165
          CALL FN15(EASP2,KE,EP,0.0,CONST)
          IF(STORE(LBNAME('TE  '))) CALL FN0(LBNAME('TE  '),EASP2)
C.... Put characteristic particle response time Tp into EASP3
          CALL FN56(EASP3,VOL,DEN2,INTFRC,1.0)
          CALL FN26(EASP3,R2)
          IF(STORE(LBNAME('TP  '))) CALL FN0(LBNAME('TP  '),EASP3)
C.... Put equivalent of 2.RHO2.R1.R2./Tp into CO
          CALL FN21(CO,R1,INTFRC,0.0,2.0)
C.... CO = GRND1 - Mostafa & Mongia
          IF(ISC.EQ.2) THEN
C.... Put Te+Tp into EASP8
            CALL FN10(EASP8,EASP2,EASP3,0.0,1.0,1.0)
C.... Now get Tp/(Te+Tp) and put in EASP3
            CALL FN27(EASP3,EASP8)
C.... Now get CO = CO * (F(EASP3+I) = 2.0*R1*R2*FIP*(Tp/(Te+Tp))
C     (I from 1 to NY*NX).
            CALL FN26(CO,EASP3)
C.... CO = GRND2 - Chen & Wood
          ELSE
            IF(INDVAR.EQ.KE) THEN
              L0CO=L0F(CO)
              L0TE=L0F(EASP2)
              L0TP=L0F(EASP3)
              DO 20 IX=IXF,IXL
                DO 20 IY=IYF,IYL
                  I=IY+NY*(IX-1)
                  F(L0CO+I)=F(L0CO+I)*(1.0-EXP(-0.0825*F(L0TP+I)
     1                                            /(F(L0TE+I)+TINY)))
  20          CONTINUE
            ENDIF
          ENDIF
C.... VAL = GRND3 - Bubble-induced turbulence production
        ELSEIF(ISC.EQ.15) THEN
          IF(INDVAR.EQ.KE) THEN
            IF(EQZ(EL1A)) EL1A=0.05
            CALL FNVSLP(1,L0F(JPRODB),CFIPA)
cc            call prn('pr 1',jprodb)
            CALL FN55(VAL,INTFRC,JPRODB,JPRODB,EL1A)
            CALL FN0(JPRODB,VAL)
            IF(STORE(LBNAME('PRKB'))) CALL FN0(LBNAME('PRKB'),JPRODB)
cc            call prn('pr 2',jprodb)
          ELSEIF(INDVAR.EQ.EP) THEN
            CALL FN56(VAL,JPRODB,EP,KE,C1E)
cc            call prn('epvl',val)
          ENDIF
          CALL FN26(VAL,R1)
        ENDIF
      ENDIF
      NAMSUB='gxdisp'
      END
c
C**** SUBROUTINE GXIPST is called from GREX3, group 13, and is entered
C     only when the patch name begins with character 'IPST'.
C
      SUBROUTINE GXIPST(INTFRC,CINT,NPATCH)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
     1       ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
      COMMON/IDATA/NX,NY,IFL2(118)
      COMMON/GENI/NXNY,IGFIL(59)
      CHARACTER*(*) NPATCH
      INTEGER VAL,CO,WALDIS,PATGEO
      LOGICAL FLUID1,CNV,IPT,NEZ,GTZ
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      NAMSUB = 'GXIPST'
C
      IF(CNV(INDVAR).AND.IPT(INDVAR) .AND. CINT.GT.0.0) THEN
C.... Obtain zero-location indices
        L0IP = L0F(INTFRC)
        L0VAL = L0F(VAL)
        L0VAR = L0F(INDVAR)
        L0VOLD = L0F(OLD(INDVAR))
        IF(FLUID1(INDVAR)) THEN
C.... Variable belongs to first phase
          INDPRT = INDVAR + 1
          L0U = L0F(U1)
        ELSE
C.... Variable belongs to second phase
          INDPRT = INDVAR - 1
          L0U = L0F(U2)
        ENDIF
        L0PRT = L0F(INDPRT)
        L0POLD = L0F(OLD(INDPRT))
        FACTOR = 0.25*CINT
        IF(NPATCH(5:6).EQ.'FE') THEN
C.... Fully explicit
          DO 10 I = 1,NXNY
            F(L0VAL+I) = 0.0
            IF(NEZ(F(L0U+I))) THEN
              INEXT = I + NY
              IF(GTZ(F(L0U+I))) INEXT = I - NY
              F(L0VAL+I) = FACTOR*F(L0IP+I)*
     1                     (-2.0* (F(L0PRT+I)-F(L0VAR+I))+F(L0POLD+I)-
     1                     F(L0VOLD+I)+F(L0POLD+INEXT)-F(L0VOLD+INEXT))
            ENDIF
   10     CONTINUE
        ELSEIF(NPATCH(5:6).EQ.'CN') THEN
C.... Crank-Nicholson
          FACTOR = 0.5*FACTOR
          DO 20 I = 1,NXNY
            F(L0VAL+I) = 0.0
            IF(NEZ(F(L0U+I))) THEN
              INEXT = I + NY
              IF(GTZ(F(L0U+I))) INEXT = I - NY
              F(L0VAL+I) = FACTOR*F(L0IP+I)*
     1                     (-3.0* (F(L0PRT+I)-F(L0VAR+I))+
     1                     F(L0PRT+INEXT)-F(L0VAR+INEXT)+F(L0POLD+I)-
     1                     F(L0VOLD+I)+F(L0POLD+INEXT)-F(L0VOLD+INEXT))
            ENDIF
   20     CONTINUE
        ELSE
          DO 30 I = 1,NXNY
            F(L0VAL+I) = 0.0
            IF(NEZ(F(L0U+I))) THEN
              INEXT = I + NY
              IF(GTZ(F(L0U+I))) INEXT = I - NY
              F(L0VAL+I) = FACTOR*F(L0IP+I)*
     1        (-F(L0PRT+I)+F(L0VAR+I)+F(L0PRT+INEXT)-F(L0VAR+INEXT))
            ENDIF
   30     CONTINUE
        ENDIF
      ENDIF
      NAMSUB = 'gxipst'
      END
c