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