1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CONV_ROP_g C 4 ! Purpose: Determine convection terms for gas continuity C 5 ! equation. Master routine. C 6 ! C 7 ! Author: M. Syamlal Date: 19-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_G(A_M, B_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 !----------------------------------------------- 36 37 IF (DISCRETIZE(1) == 0) THEN ! 0 & 1 => first order upwinding 38 CALL CONV_ROP_G0 (A_M) 39 ELSE 40 CALL CONV_ROP_G1 (A_M) 41 ENDIF 42 43 44 RETURN 45 END SUBROUTINE CONV_ROP_G 46 47 48 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 49 ! C 50 ! Subroutine: CONV_ROP_g0 C 51 ! Purpose: Determine convection terms for gas continuity equation C 52 ! C 53 ! Notes: The off-diagonal coefficients calculated here must be C 54 ! positive. The center coefficient and the source vector C 55 ! are negative. -- First order upwinding C 56 ! C 57 ! Author: M. Syamlal Date: 2-JUL-96 C 58 ! Reviewer: Date: C 59 ! C 60 ! Literature/Document References: C 61 ! Variables referenced: C 62 ! Variables modified: C 63 ! Local variables: C 64 ! C 65 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 66 67 SUBROUTINE CONV_ROP_G0(A_M) 68 69 !----------------------------------------------- 70 ! Modules 71 !----------------------------------------------- 72 USE param 73 USE param1 74 USE fldvar 75 USE run 76 USE parallel 77 USE physprop 78 USE geometry 79 USE indices 80 USE pgcor 81 USE compar 82 USE functions 83 IMPLICIT NONE 84 !----------------------------------------------- 85 ! Dummy arguments 86 !----------------------------------------------- 87 ! Septadiagonal matrix A_m 88 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 89 !----------------------------------------------- 90 ! Local variables 91 !----------------------------------------------- 92 ! Indices 93 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 94 INTEGER :: IMJK, IJMK, IJKM 95 !----------------------------------------------- 96 97 98 ! Calculate convection fluxes through each of the faces 99 100 !!$omp parallel do private( J, K, IJK, IPJK, IJPK, IJKP, & 101 !!$omp& IMJK, IJMK, IJKM) & 102 !!$omp& schedule(static) 103 DO IJK = ijkstart3, ijkend3 104 IF (PHASE_4_P_G(IJK) /= 0) THEN 105 ! if phase_4_p_g is assigned the index of a solids phase then this 106 ! section is entered. currently, phase_4_p_g is always assigned to 107 ! the gas phase index 0 and so this section will never be entered. 108 109 ! it was previously possible for phase_4_p_g to be assigned the index of 110 ! the gas phase in some cells, and the index of a solids phase in other 111 ! cells if that solids phase has close_packed=F and was in higher 112 ! concentrations than the gas. in such a case this branch would become 113 ! activated, while the corresponding solids continuity would become 114 ! skipped. moreover, that solids phase's volume fraction would be 115 ! corrected rather than the gas phases void fraction in calc_vol_fr. 116 I = I_OF(IJK) 117 J = J_OF(IJK) 118 K = K_OF(IJK) 119 IPJK = IP_OF(IJK) 120 IJPK = JP_OF(IJK) 121 IJKP = KP_OF(IJK) 122 IMJK = IM_OF(IJK) 123 IJMK = JM_OF(IJK) 124 IJKM = KM_OF(IJK) 125 126 ! East face (i+1/2, j, k) 127 A_M(IJK,east,0) = ZMAX((-U_G(IJK)))*AYZ(IJK) 128 A_M(IPJK,west,0) = ZMAX(U_G(IJK))*AYZ(IJK) 129 130 ! North face (i, j+1/2, k) 131 A_M(IJK,north,0) = ZMAX((-V_G(IJK)))*AXZ(IJK) 132 A_M(IJPK,south,0) = ZMAX(V_G(IJK))*AXZ(IJK) 133 134 ! Top face (i, j, k+1/2) 135 IF (DO_K) THEN 136 A_M(IJK,top,0) = ZMAX((-W_G(IJK)))*AXY(IJK) 137 A_M(IJKP,bottom,0) = ZMAX(W_G(IJK))*AXY(IJK) 138 ENDIF 139 140 ! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2) 141 ! faces if the neighboring west, south, bottom cells have 142 ! phase_4_p_g of 0 143 IF(PHASE_4_P_G(IMJK)==0) A_M(IJK,west,0)=ZMAX(U_G(IMJK))*AYZ(IMJK) 144 IF(PHASE_4_P_G(IJMK)==0) A_M(IJK,south,0)=ZMAX(V_G(IJMK))*AXZ(IJMK) 145 IF (DO_K) THEN 146 IF(PHASE_4_P_G(IJKM)==0) A_M(IJK,bottom,0)=ZMAX(W_G(IJKM))*AXY(IJKM) 147 ENDIF 148 149 ENDIF 150 ENDDO 151 152 RETURN 153 END SUBROUTINE CONV_ROP_G0 154 155 156 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 157 ! C 158 ! Subroutine: CONV_ROP_g1 C 159 ! Purpose: Determine convection terms for gas continuity equation. C 160 ! C 161 ! Notes: The off-diagonal coefficients calculated here must be C 162 ! positive. The center coefficient and the source vector C 163 ! are negative. -- Higher order schemes C 164 ! C 165 ! Author: M. Syamlal Date:19-MAR-97 C 166 ! Reviewer: Date: C 167 ! C 168 ! Literature/Document References: C 169 ! Variables referenced: C 170 ! Variables modified: C 171 ! Local variables: C 172 ! C 173 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 174 175 SUBROUTINE CONV_ROP_G1(A_M) 176 177 !----------------------------------------------- 178 ! Modules 179 !----------------------------------------------- 180 USE compar 181 USE fldvar 182 USE functions 183 USE geometry 184 USE indices 185 USE parallel 186 USE param 187 USE param1 188 USE pgcor 189 USE physprop 190 USE run 191 USE xsi 192 IMPLICIT NONE 193 !----------------------------------------------- 194 ! D u m m y A r g u m e n t s 195 !----------------------------------------------- 196 ! Septadiagonal matrix A_m 197 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 198 !----------------------------------------------- 199 ! Local variables 200 !----------------------------------------------- 201 ! Indices 202 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP 203 INTEGER :: IMJK, IJMK, IJKM 204 ! loezos 205 INTEGER :: incr 206 ! temporary use of global arrays: 207 ! xsi_array: convection weighting factors 208 DOUBLE PRECISION :: XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),& 209 XSI_t(DIMENSION_3) 210 !----------------------------------------------- 211 212 213 ! loezos 214 incr=0 215 216 ! Calculate convection factors 217 CALL CALC_XSI (DISCRETIZE(1), ROP_G, U_G, V_G, W_G, XSI_E, & 218 XSI_N, XSI_T,incr) 219 220 ! Calculate convection fluxes through each of the faces 221 222 !!$omp parallel do private( J, K, IJK, IPJK, IJPK, IJKP, & 223 !!$omp& IMJK, IJMK, IJKM) & 224 !!$omp& schedule(static) 225 DO IJK = ijkstart3, ijkend3 226 IF (PHASE_4_P_G(IJK) /= 0) THEN 227 I = I_OF(IJK) 228 J = J_OF(IJK) 229 K = K_OF(IJK) 230 IPJK = IP_OF(IJK) 231 IJPK = JP_OF(IJK) 232 IJKP = KP_OF(IJK) 233 IMJK = IM_OF(IJK) 234 IJMK = JM_OF(IJK) 235 IJKM = KM_OF(IJK) 236 237 ! East face (i+1/2, j, k) 238 A_M(IJK,east,0) = -XSI_E(IJK)*U_G(IJK)*AYZ(IJK) 239 A_M(IPJK,west,0) = (ONE - XSI_E(IJK))*U_G(IJK)*AYZ(IJK) 240 241 ! North face (i, j+1/2, k) 242 A_M(IJK,north,0) = -XSI_N(IJK)*V_G(IJK)*AXZ(IJK) 243 A_M(IJPK,south,0) = (ONE - XSI_N(IJK))*V_G(IJK)*AXZ(IJK) 244 245 ! Top face (i, j, k+1/2) 246 IF (DO_K) THEN 247 A_M(IJK,top,0) = -XSI_T(IJK)*W_G(IJK)*AXY(IJK) 248 A_M(IJKP,bottom,0) = (ONE - XSI_T(IJK))*W_G(IJK)*AXY(IJK) 249 ENDIF 250 251 ! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2) 252 ! faces if the neighboring west, south, bottom cells have 253 ! phase_4_p_g of 0 254 IF (PHASE_4_P_G(IMJK) == 0) A_M(IJK,west,0) = (ONE - XSI_E(IMJK))*U_G(& 255 IMJK)*AYZ(IMJK) 256 IF (PHASE_4_P_G(IJMK) == 0) A_M(IJK,south,0) = (ONE - XSI_N(IJMK))*V_G(& 257 IJMK)*AXZ(IJMK) 258 IF (DO_K) THEN 259 IF (PHASE_4_P_G(IJKM) == 0) A_M(IJK,bottom,0) = (ONE - XSI_T(IJKM))*& 260 W_G(IJKM)*AXY(IJKM) 261 ENDIF 262 263 ENDIF 264 ENDDO 265 266 RETURN 267 END SUBROUTINE CONV_ROP_G1 268 269 270