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