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