1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CONV_ROP_s C 4 ! Purpose: Determine convection terms for solids continuity C 5 ! equation. Master routine. C 6 ! C 7 ! Author: M. Syamlal Date: 18-MAR-97 C 8 ! Reviewer: Date: C 9 ! C 10 ! Literature/Document References: C 11 ! Variables referenced: C 12 ! Variables modified: C 13 ! Local variables: C 14 ! C 15 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 16 17 SUBROUTINE CONV_ROP_S(A_M, B_M, M) 18 19 !----------------------------------------------- 20 ! Modules 21 !----------------------------------------------- 22 USE param 23 USE param1 24 USE fldvar 25 USE run 26 USE compar 27 IMPLICIT NONE 28 !----------------------------------------------- 29 ! Dummy arguments 30 !----------------------------------------------- 31 ! Septadiagonal matrix A_m 32 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 33 ! Vector b_m 34 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M) 35 ! Solids phase index 36 INTEGER, INTENT(IN) :: M 37 !----------------------------------------------- 38 39 IF (DISCRETIZE(2) == 0) THEN ! 0 & 1 => first order upwinding 40 CALL CONV_ROP_S0 (A_M, M) 41 ELSE 42 CALL CONV_ROP_S1 (A_M, M) 43 ENDIF 44 45 RETURN 46 END SUBROUTINE CONV_ROP_S 47 48 49 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 50 ! C 51 ! Subroutine: CONV_ROP_s0 C 52 ! Purpose: Determine convection terms for solids continuity equation. C 53 ! C 54 ! Notes: The off-diagonal coefficients calculated here must be C 55 ! positive. The center coefficient and the source vector C 56 ! are negative. -- First order upwinding C 57 ! C 58 ! Author: M. Syamlal Date: 3-JUL-96 C 59 ! Reviewer: Date: C 60 ! C 61 ! Literature/Document References: C 62 ! Variables referenced: C 63 ! Variables modified: C 64 ! Local variables: C 65 ! C 66 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 67 68 SUBROUTINE CONV_ROP_S0(A_M, M) 69 70 !----------------------------------------------- 71 ! Modules 72 !----------------------------------------------- 73 USE param 74 USE param1 75 USE fldvar 76 USE run 77 USE parallel 78 USE physprop 79 USE geometry 80 USE indices 81 USE pgcor 82 USE pscor 83 USE compar 84 USE functions 85 IMPLICIT NONE 86 !----------------------------------------------- 87 ! Dummy arguments 88 !----------------------------------------------- 89 ! Septadiagonal matrix A_m 90 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 91 ! Solids phase index 92 INTEGER, INTENT(IN) :: M 93 !----------------------------------------------- 94 ! Local variables 95 !----------------------------------------------- 96 ! Indices 97 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 98 INTEGER :: IMJK, IJMK, IJKM 99 !----------------------------------------------- 100 101 ! Calculate convection-diffusion fluxes through each of the faces 102 103 !!$omp parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP, & 104 !!$omp& IMJK, IJMK, IJKM) & 105 !!$omp& schedule(static) 106 DO IJK = ijkstart3, ijkend3 107 IF (PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN 108 ! if either phase_4_p_g or phase_4_p_s are assigned the index of the 109 ! current solids phase then this section is skipped. currently, 110 ! phase_4_p_g and phase_4_p_s are always assigned to 0 and undefined, 111 ! respectively. Hence this section should always be entered. 112 113 ! it was previously possible for phase_4_p_g to be assigned the index of 114 ! the gas phase in some cells, and the index of a solids phase in other 115 ! cells if that solids phase has close_packed=F and was in higher 116 ! concentrations than the gas. in such a case this branch would become 117 ! skipped, while the corresponding gas continuity would become 118 ! activated. moreover, that solids phase's volume fraction would be 119 ! corrected rather than the gas phases void fraction in calc_vol_fr. 120 121 ! if phase_4_p_s were to be assigned the index of a solids phase that 122 ! can close pack if the current cell exhibited close packing, then 123 ! this branch would become skipped. moreover, that solids phase's 124 ! volume fraction would then be corected based on the value of 125 ! maximum packing in that cell and the sum of all other solids that 126 ! can close pack in calc_vol_fr. 127 I = I_OF(IJK) 128 J = J_OF(IJK) 129 K = K_OF(IJK) 130 IPJK = IP_OF(IJK) 131 IJPK = JP_OF(IJK) 132 IJKP = KP_OF(IJK) 133 IMJK = IM_OF(IJK) 134 IJMK = JM_OF(IJK) 135 IJKM = KM_OF(IJK) 136 137 ! East face (i+1/2, j, k) 138 A_M(IJK,east,M) = ZMAX((-U_S(IJK,M)))*AYZ(IJK) 139 A_M(IPJK,west,M) = ZMAX(U_S(IJK,M))*AYZ(IJK) 140 141 ! North face (i, j+1/2, k) 142 A_M(IJK,north,M) = ZMAX((-V_S(IJK,M)))*AXZ(IJK) 143 A_M(IJPK,south,M) = ZMAX(V_S(IJK,M))*AXZ(IJK) 144 145 ! Top face (i, j, k+1/2) 146 IF (DO_K) THEN 147 A_M(IJK,top,M) = ZMAX((-W_S(IJK,M)))*AXY(IJK) 148 A_M(IJKP,bottom,M) = ZMAX(W_S(IJK,M))*AXY(IJK) 149 ENDIF 150 151 ! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2) 152 ! faces if the neighboring west, south, bottom cells have 153 ! phase_4_p_g of m or phase_4_p_s of m. 154 IF (PHASE_4_P_G(IMJK)==M .OR. PHASE_4_P_S(IMJK)==M) & 155 A_M(IJK,west,M) = ZMAX(U_S(IMJK,M))*AYZ(IMJK) 156 IF (PHASE_4_P_G(IJMK)==M .OR. PHASE_4_P_S(IJMK)==M) & 157 A_M(IJK,south,M) = ZMAX(V_S(IJMK,M))*AXZ(IJMK) 158 159 IF (DO_K) THEN 160 IF (PHASE_4_P_G(IJKM)==M .OR. PHASE_4_P_S(IJKM)==M) & 161 A_M(IJK,bottom,M) = ZMAX(W_S(IJKM,M))*AXY(IJKM) 162 ENDIF 163 ENDIF 164 END DO 165 166 RETURN 167 END SUBROUTINE CONV_ROP_S0 168 169 170 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 171 ! C 172 ! Subroutine: CONV_ROP_s1 C 173 ! Purpose: Determine convection terms for solids continuity equation. C 174 ! Notes: The off-diagonal coefficients calculated here must be C 175 ! positive. The center coefficient and the source vector C 176 ! are negative. -- Higher order methods C 177 ! See conv_rop_s0 for additional details. C 178 ! C 179 ! Author: M. Syamlal Date: 18-MAR-97 C 180 ! Reviewer: Date: C 181 ! C 182 ! C 183 ! Literature/Document References: C 184 ! C 185 ! Variables referenced: C 186 ! Variables modified: C 187 ! C 188 ! Local variables: C 189 ! C 190 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 191 192 SUBROUTINE CONV_ROP_S1(A_M, M) 193 194 !----------------------------------------------- 195 ! Modules 196 !----------------------------------------------- 197 USE compar 198 USE fldvar 199 USE functions 200 USE geometry 201 USE indices 202 USE parallel 203 USE param 204 USE param1 205 USE pgcor 206 USE physprop 207 USE pscor 208 USE run 209 USE xsi 210 IMPLICIT NONE 211 !----------------------------------------------- 212 ! Dummy arguments 213 !----------------------------------------------- 214 ! Septadiagonal matrix A_m 215 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 216 ! Solids phase index 217 INTEGER, INTENT(IN) :: M 218 !----------------------------------------------- 219 ! Local variables 220 !----------------------------------------------- 221 ! Indices 222 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 223 INTEGER :: IMJK, IJMK, IJKM 224 ! loezos 225 INTEGER :: incr 226 ! temporary use of global arrays: 227 ! xsi_array: convection weighting factors 228 DOUBLE PRECISION :: XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),& 229 XSI_t(DIMENSION_3) 230 !----------------------------------------------- 231 232 ! Loezos: 233 incr=0 234 235 ! Calculate convection factors 236 CALL CALC_XSI(DISCRETIZE(2), ROP_S(1,M), U_S(1,M), V_S(1,M),& 237 W_S(1,M), XSI_E, XSI_N, XSI_T, incr) 238 239 ! Calculate convection-diffusion fluxes through each of the faces 240 241 !!$omp parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP, & 242 !!$omp& IMJK, IJMK, IJKM) & 243 !!$omp& schedule(static) 244 DO IJK = ijkstart3, ijkend3 245 IF (PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN 246 I = I_OF(IJK) 247 J = J_OF(IJK) 248 K = K_OF(IJK) 249 IPJK = IP_OF(IJK) 250 IJPK = JP_OF(IJK) 251 IJKP = KP_OF(IJK) 252 IMJK = IM_OF(IJK) 253 IJMK = JM_OF(IJK) 254 IJKM = KM_OF(IJK) 255 256 ! East face (i+1/2, j, k) 257 A_M(IJK,east,M) = -XSI_E(IJK)*U_S(IJK,M)*AYZ(IJK) 258 A_M(IPJK,west,M) = (ONE - XSI_E(IJK))*U_S(IJK,M)*AYZ(IJK) 259 260 ! North face (i, j+1/2, k) 261 A_M(IJK,north,M) = -XSI_N(IJK)*V_S(IJK,M)*AXZ(IJK) 262 A_M(IJPK,south,M) = (ONE - XSI_N(IJK))*V_S(IJK,M)*AXZ(IJK) 263 264 ! Top face (i, j, k+1/2) 265 IF (DO_K) THEN 266 A_M(IJK,top,M) = -XSI_T(IJK)*W_S(IJK,M)*AXY(IJK) 267 A_M(IJKP,bottom,M) = (ONE - XSI_T(IJK))*W_S(IJK,M)*AXY(IJK) 268 ENDIF 269 270 IF (PHASE_4_P_G(IMJK)==M .OR. PHASE_4_P_S(IMJK)==M) A_M(IJK,west,M) = & 271 (ONE - XSI_E(IMJK))*U_S(IMJK,M)*AYZ(IMJK) 272 IF (PHASE_4_P_G(IJMK)==M .OR. PHASE_4_P_S(IJMK)==M) A_M(IJK,south,M) = & 273 (ONE - XSI_N(IJMK))*V_S(IJMK,M)*AXZ(IJMK) 274 IF (DO_K) THEN 275 IF (PHASE_4_P_G(IJKM)==M .OR. PHASE_4_P_S(IJKM)==M) A_M(IJK,bottom,M)& 276 = (ONE - XSI_T(IJKM))*W_S(IJKM,M)*AXY(IJKM) 277 ENDIF 278 ENDIF 279 280 ENDDO 281 282 RETURN 283 END SUBROUTINE CONV_ROP_S1 284 285