1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CONV_Pp_g C 4 ! Purpose: Determine convection terms for Pressure correction C 5 ! equation. C 6 ! C 7 ! Notes: The off-diagonal coefficients calculated here must be C 8 ! positive. The center coefficient and the source vector are C 9 ! negative. C 10 ! Multiplication with factors d_e, d_n, and d_t are carried C 11 ! out in source_pp_g. Constant pressure boundaries are C 12 ! handled by holding the Pp_g at the boundaries zero. For C 13 ! specified mass flow boundaries (part of) a's are calculated C 14 ! here since b is calculated from a's in source_pp_g. After C 15 ! calculating b, a's are multiplied by d and at the flow C 16 ! boundaries are set to zero. C 17 ! C 18 ! Author: M. Syamlal Date: 20-JUN-96 C 19 ! Reviewer: Date: C 20 ! C 21 ! Revision Number: 1 C 22 ! Purpose: to use face densities calculated in CONV_ROP C 23 ! Author: M. Syamlal Date: 1-JUN-05 C 24 ! Reviewer: Date: C 25 ! C 26 ! Literature/Document References: C 27 ! Variables referenced: C 28 ! Variables modified: C 29 ! Local variables: C 30 ! C 31 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 32 33 SUBROUTINE CONV_PP_G(A_M, B_M) 34 35 !----------------------------------------------- 36 ! Modules 37 !----------------------------------------------- 38 USE param 39 USE param1 40 USE fldvar 41 USE run 42 USE parallel 43 USE physprop 44 USE geometry 45 USE indices 46 USE pgcor 47 USE compar 48 USE mflux 49 USE functions 50 IMPLICIT NONE 51 !----------------------------------------------- 52 ! Dummy arguments 53 !----------------------------------------------- 54 ! Septadiagonal matrix A_m 55 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 56 ! Vector b_m 57 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M) 58 !----------------------------------------------- 59 ! Local variables 60 !----------------------------------------------- 61 ! Indices 62 INTEGER :: IJK, IPJK, IJPK, IJKP 63 INTEGER :: IMJK, IJMK, IJKM 64 INTEGER :: M 65 ! local value of A_m 66 DOUBLE PRECISION :: am 67 !----------------------------------------------- 68 69 70 ! Calculate convection fluxes through each of the faces 71 72 !$omp parallel do default(none) private( IJK, IPJK, IJPK, IJKP, M, AM, IMJK, IJMK, IJKM) & 73 !$omp shared(ijkstart3, ijkend3, rop_ge, rop_gn, rop_gt, axz, axy, ayz, do_k, a_m) 74 DO IJK = ijkstart3, ijkend3 75 IF (FLUID_AT(IJK)) THEN 76 IPJK = IP_OF(IJK) 77 IJPK = JP_OF(IJK) 78 IJKP = KP_OF(IJK) 79 80 ! East face (i+1/2, j, k) 81 AM = ROP_GE(IJK)*AYZ(IJK) 82 A_M(IJK,east,0) = AM 83 A_M(IPJK,west,0) = AM 84 85 ! North face (i, j+1/2, k) 86 AM = ROP_GN(IJK)*AXZ(IJK) 87 A_M(IJK,north,0) = AM 88 A_M(IJPK,south,0) = AM 89 90 ! Top face (i, j, k+1/2) 91 IF (DO_K) THEN 92 AM = ROP_GT(IJK)*AXY(IJK) 93 A_M(IJK,top,0) = AM 94 A_M(IJKP,bottom,0) = AM 95 ENDIF 96 97 ! West face (i-1/2, j, k) 98 IMJK = IM_OF(IJK) 99 IF (.NOT.FLUID_AT(IMJK)) THEN 100 AM = ROP_GE(IMJK)*AYZ(IMJK) 101 A_M(IJK,west,0) = AM 102 ENDIF 103 104 ! South face (i, j-1/2, k) 105 IJMK = JM_OF(IJK) 106 IF (.NOT.FLUID_AT(IJMK)) THEN 107 AM = ROP_GN(IJMK)*AXZ(IJMK) 108 A_M(IJK,south,0) = AM 109 ENDIF 110 111 ! Bottom face (i, j, k-1/2) 112 IF (DO_K) THEN 113 IJKM = KM_OF(IJK) 114 IF (.NOT.FLUID_AT(IJKM)) THEN 115 AM = ROP_GT(IJKM)*AXY(IJKM) 116 A_M(IJK,bottom,0) = AM 117 ENDIF 118 ENDIF 119 ENDIF ! end if (fluid_at(ijk)) 120 ENDDO ! end do ijk=ijkstart3,ijkend3 121 122 123 ! solids phase terms are needed for those solids phases that do not 124 ! become close packed. these terms are used to essentially make the 125 ! matrix equation for gas pressure correction a mixture pressure 126 ! correction equation by adding solids phases that do not become 127 ! close packed 128 DO M = 1, MMAX 129 IF (.NOT.CLOSE_PACKED(M)) THEN 130 DO IJK = ijkstart3, ijkend3 131 IF (FLUID_AT(IJK)) THEN 132 IPJK = IP_OF(IJK) 133 IJPK = JP_OF(IJK) 134 IJKP = KP_OF(IJK) 135 136 ! East face (i+1/2, j, k) 137 AM = ROP_SE(IJK,M)*AYZ(IJK) 138 A_M(IJK,east,M) = AM 139 A_M(IPJK,west,M) = AM 140 141 ! North face (i, j+1/2, k) 142 AM = ROP_SN(IJK,M)*AXZ(IJK) 143 A_M(IJK,north,M) = AM 144 A_M(IJPK,south,M) = AM 145 146 ! Top face (i, j, k+1/2) 147 IF (DO_K) THEN 148 AM = ROP_ST(IJK,M)*AXY(IJK) 149 A_M(IJK,top,M) = AM 150 A_M(IJKP,bottom,M) = AM 151 ENDIF 152 153 ! West face (i-1/2, j, k) 154 IMJK = IM_OF(IJK) 155 IF (.NOT.FLUID_AT(IMJK)) THEN 156 AM = ROP_SE(IMJK,M)*AYZ(IMJK) 157 A_M(IJK,west,M) = AM 158 ENDIF 159 160 ! South face (i, j-1/2, k) 161 IJMK = JM_OF(IJK) 162 IF (.NOT.FLUID_AT(IJMK)) THEN 163 AM = ROP_SN(IJMK,M)*AXZ(IJMK) 164 A_M(IJK,south,M) = AM 165 ENDIF 166 167 ! Bottom face (i, j, k-1/2) 168 IF (DO_K) THEN 169 IJKM = KM_OF(IJK) 170 IF (.NOT.FLUID_AT(IJKM)) THEN 171 AM = ROP_ST(IJKM,M)*AXY(IJKM) 172 A_M(IJK,bottom,M) = AM 173 ENDIF 174 ENDIF 175 176 ENDIF ! end if(fluid_at(ijk)) 177 ENDDO ! end do ijk=ijkstart3,ijkend3 178 ENDIF ! end if(.not.close_packed) 179 ENDDO ! end do (m=1,mmax) 180 181 RETURN 182 END SUBROUTINE CONV_PP_G 183