cgxidif.htm c c
c... File name ... GXIDIF.HTM ... 050713
C#### JCL 05.07.13 remove call to PRPADR, change .eq.to == etc.
C#### JCL 08.09.10 initialise INTDIF in case not set by function, otherwise
C####              undefined e.g. when CINT()=GRND10 for InForm
      REAL FUNCTION INTDIF(I)
      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/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(39),
     1            ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IGF4(4)
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      LOGICAL GRN
c-----------------------------------------------------------------------
      INTDIF=-999.0
C.... Bulk-to-interface diffusive-transfer coefficient:
      IF(IGRND==-1) THEN
        INTDIF=PRPRTY
      ELSEIF(IGRND==1) THEN
        INTDIF=0.0
      ELSEIF(IGRND==7 .OR. IGRND==8) THEN
C.... Use stored slip velocity, or calculate it (CINH1C stores minimum
C     slip velocity):
        IF(STVSLP) THEN
          VSLP= F(L0VSLP+I)
        ELSE
          VSLP= AMAX1(CINH1C, VSLPCL(I,XCYCZ,SLUNX1))
          IF(LBVREL/=0) F(L0VSLP+I)= VSLP
        ENDIF
        IF(L0DIAM/=0.AND.L0REYD/=0) THEN
          DIAM= F(L0DIAM+I); REYD= F(L0REYD+I)
        ELSE
C.... Compute particle size Dp= Dp0*(Rd/RS)**0.3333, where Dp0=CINH2C
C     is the initial particle diameter for the SHADOW method, and
C     particle Reynolds number ReD= Dp*RelVel/Enul:
          DIAM= ABS(CINH2C)
          IF(SOLVRS)
     1       DIAM= DIAM*AMAX1(0., AMIN1(1.,F(L0RD+I)/F(L0RS+I)))**0.3333
          REYD= DIAM*VSLP/(F(L0VISL+I)+TINY)
          F(L0DIAM+I)= DIAM; F(L0REYD+I)= REYD
        ENDIF
C.... Get Prandtl number and diffusivity:
        PRLM= PRNDTL(INDVAR)
        IF(GRN(-ABS(PRLM))) THEN
C.... GROUND-set Prandtl number (PRLM=GRND) or diffusivity (PRLM=-GRND):
          IF(PRLM<0.0) THEN
            PRLM= F(L0PRL+I); DIFF= F(L0VISL+I)/PRLM
          ELSE
            DIFF= F(L0PRL+I); PRLM= F(L0VISL+I)/DIFF
          ENDIF
        ELSE
C.... Constant Prandtl (PRLM>0.0) or diffusivity (PRLM<0.0):
          IF(PRLM>0.0) THEN
            DIFF= F(L0VISL+I)/PRLM
          ELSE
            DIFF= - PRLM; PRLM= F(L0VISL+I)/DIFF
          ENDIF
        ENDIF
        IF(IGRND==7) THEN
C.... Calculate Nusselt number from Clift-Grace-Weber correlation:
          FAC= ((1.0 + 1.0/(REYD*PRLM))*PRLM)**0.3333
          IF(REYD<=400) THEN
            DNUSS= 1.0 + FAC*REYD**0.41
          ELSEIF(REYD<=2000.) THEN
            DNUSS= 1.0 + 0.752*FAC*REYD**0.472
          ELSE
            DNUSS= 1.0 + FAC*(0.44*SQRT(REYD) + 0.034*REYD**0.71)
          ENDIF
        ELSE
C.... Calculate Nusselt number from Ranz-Marshall correlation:
          DNUSS= 2.0 + 0.6*SQRT(REYD)*PRLM**0.3333
        ENDIF
        IF(LBNUSS/=0) F(L0NUSS+I)= DNUSS
C.... Compute the surface area of particles (spheres) in each cell:
        ASD= 6.*F(L0RD+I)*CELVOL(I)/DIAM
C.... Set interface transfer coefficient Ci= RHOCG*(Enul*Nu/Prl)/Dp*As:
        INTDIF= CINHFG*F(L0DEN+I)*DIFF*DNUSS/DIAM*ASD
      ELSEIF(IGRND==9) THEN
C.... Constant diffusive-transfer coefficient:
        INTDIF= CINHAG
      ENDIF
      END
C-----------------------------------------------------------------------
      SUBROUTINE SLBIDF(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
C#### SCM 04.02.09 Split commons for REALs and INTEGERs
C     1      /TSKEM/ GKTDKP,LBKP,LBKT,LBET,LBVOSQ,LBOMEG
      COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
     1      /TSKEMR/ GKTDKP
     1      /VELCMN/L0UVW(6),L0UVW2(6)
     1      /CELPAR/IPHASE,IPROP,IGRND,IFILEP,KPROP
      COMMON/GENI/IGF1(2),NXNYST,NDIR,KDUMM,IGF2(4),NFM,IGF3(21),IPRL,
     1           IBTAU,IGF4(16),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
     1           IPRPS,IRADX,IRADY,IRADZ,IVFOL
C#### JCL 05.07.13 add KZXCY
     1      /F0/ IF01(29),L0XYDX,L0XYDY,IF02(3),L0XYRV,L0XYXG,IF03,
     1           L0XYYG,IF04,L0XYDXG,L0XYDYG,IF05(68),KZXCY,IF05A(37),
     1           L0AHZ,IF06(17),L0XYDZ,L0XYDZG,IF07(137)
      COMMON/NAMFN/NAMFUN,NAMSUB
      REAL INTDIF
      LOGICAL DBGLOC,SLD,GRN,CINT78,FL1CON
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB= 'SLBIDF'
      if(flag.or.dbgloc) call banner(1,namsub,050713)
      IGR= 10
      ISC= IPROP-20
C.... Call GROUND for the user set property:
      IF(IGRND==0) THEN
        IF(USEGRD) THEN
          CALL GROUND
        ENDIF
        GO TO 800
      ENDIF
C.... For a constant property go to the slab loop
      IF(IGRND==-1) GO TO 700
C#### JCL 05.07.13 add needed settings here      
C####      CALL PRPADR
C.... Set constants and other auxiliary variables:
C-----------------------------------------------------------------------
C.... Bulk-to-interface diffusive-transfer coeff.:
        IF(IPHASE==1) THEN
          CINHAG = CINH1A; CINHAG = CINH2A
        ENDIF
        CINT78= .FALSE.
        IF(IPILOPT==7 .OR. IPILOPT==8) THEN
          L0VISL= L0F(VISL)
          IF(GRN(-ABS(PRNDTL(INDVAR)))) L0PRL= L0F(IPRL)
          IF(LBSIZE/=0) L0DIAM= L0F(LBSIZE)
          IF(LBREYD/=0) L0REYD= L0F(LBREYD)
          IF(LBNUSS/=0) L0NUSS= L0F(LBNUSS)
C.... First check the option used to calculate Cfip. If it is GRND7 or
C     GRND8, use values of slip velocity, particle sizes and Reynolds
C     number calculated for Cfip.
          IF(CFIPS>=0) THEN
            IGRNDCF=-1
          ELSE  
            IGRNDCF= NINT(ABS(CFIPS-GRND))/10
          ENDIF
          CFIP78= IGRNDCF==7 .OR. IGRNDCF==8
          STVSLP= CFIP78 .AND. LBVREL/=0
          IF(STVSLP) L0VSLP= L0F(LBVREL)
C.... Determine indices for continuous and dispersed phases:
          IF(CFIP78) THEN
            IF(IGRNDCF==7) THEN
              L0DEN= L0F(DEN1); L0RD = L0R2
            ELSE
              L0DEN= L0F(DEN2); L0RD = L0R1
            ENDIF
          ELSEIF(IPHASE==1) THEN
            L0DEN= L0F(DEN1); L0RD = L0R2
          ELSE
            L0DEN= L0F(DEN2); L0RD = L0R1
          ENDIF
          INEIGH  = ITWO(INDVAR+1,INDVAR-1,IPHASE==1)
          PHNHFG= 1.0
          IF(CINT(INEIGH)>1.E10) THEN
            PHNHFG = 0.5; CINT78  = IPHASE==1
          ENDIF
          IF(XCYCLE) XCYCZ=NEZ(F(KZXCY+IZ))
          IF(SOLVE(3))  L0UVW(3)=L0F(3);  IF(SOLVE(4))  L0UVW2(3)=L0F(4)
          IF(SOLVE(5))  L0UVW(1)=L0F(5);  IF(SOLVE(6))  L0UVW2(1)=L0F(6)
          IF(SOLVE(7))  L0UVW(5)=L0F(7);  IF(SOLVE(8))  L0UVW2(5)=L0F(8)
          IF(SOLVRS) L0RS = L0F(11)
        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) = INTDIF(I)
      ELSE
C.... exclude solids
        DO 70 I= 1,NXNYST
          IF(SLD(I)) THEN
            F(KPROP+I)= TINY
          ELSE
            F(KPROP+I)= INTDIF(I)
          ENDIF
  70    CONTINUE
      ENDIF
C----------------------------------------------------------------------
C.... Corrections, debug print-out, and other property adjustments:
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
      NAMSUB= 'slbidf'
      if(flag.or.dbgloc) call banner(2,namsub,0)
      END
c