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 matrix 44 USE physprop 45 USE geometry 46 USE indices 47 USE pgcor 48 USE compar 49 USE mflux 50 USE functions 51 IMPLICIT NONE 52 !----------------------------------------------- 53 ! Dummy arguments 54 !----------------------------------------------- 55 ! Septadiagonal matrix A_m 56 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 57 ! Vector b_m 58 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M) 59 !----------------------------------------------- 60 ! Local variables 61 !----------------------------------------------- 62 ! Indices 63 INTEGER :: IJK, IPJK, IJPK, IJKP 64 INTEGER :: IMJK, IJMK, IJKM 65 INTEGER :: M 66 ! local value of A_m 67 DOUBLE PRECISION :: am 68 !----------------------------------------------- 69 70 71 ! Calculate convection fluxes through each of the faces 72 73 !$omp parallel do default(none) private( IJK, IPJK, IJPK, IJKP, M, AM, IMJK, IJMK, IJKM) & 74 !$omp shared(ijkstart3, ijkend3, rop_ge, rop_gn, rop_gt, axz, axy, ayz, do_k, a_m) 75 DO IJK = ijkstart3, ijkend3 76 IF (FLUID_AT(IJK)) THEN 77 IPJK = IP_OF(IJK) 78 IJPK = JP_OF(IJK) 79 IJKP = KP_OF(IJK) 80 81 ! East face (i+1/2, j, k) 82 AM = ROP_GE(IJK)*AYZ(IJK) 83 A_M(IJK,E,0) = AM 84 A_M(IPJK,W,0) = AM 85 86 ! North face (i, j+1/2, k) 87 AM = ROP_GN(IJK)*AXZ(IJK) 88 A_M(IJK,N,0) = AM 89 A_M(IJPK,S,0) = AM 90 91 ! Top face (i, j, k+1/2) 92 IF (DO_K) THEN 93 AM = ROP_GT(IJK)*AXY(IJK) 94 A_M(IJK,T,0) = AM 95 A_M(IJKP,B,0) = AM 96 ENDIF 97 98 ! West face (i-1/2, j, k) 99 IMJK = IM_OF(IJK) 100 IF (.NOT.FLUID_AT(IMJK)) THEN 101 AM = ROP_GE(IMJK)*AYZ(IMJK) 102 A_M(IJK,W,0) = AM 103 ENDIF 104 105 ! South face (i, j-1/2, k) 106 IJMK = JM_OF(IJK) 107 IF (.NOT.FLUID_AT(IJMK)) THEN 108 AM = ROP_GN(IJMK)*AXZ(IJMK) 109 A_M(IJK,S,0) = AM 110 ENDIF 111 112 ! Bottom face (i, j, k-1/2) 113 IF (DO_K) THEN 114 IJKM = KM_OF(IJK) 115 IF (.NOT.FLUID_AT(IJKM)) THEN 116 AM = ROP_GT(IJKM)*AXY(IJKM) 117 A_M(IJK,B,0) = AM 118 ENDIF 119 ENDIF 120 ENDIF ! end if (fluid_at(ijk)) 121 ENDDO ! end do ijk=ijkstart3,ijkend3 122 123 124 ! solids phase terms are needed for those solids phases that do not 125 ! become close packed. these terms are used to essentially make the 126 ! matrix equation for gas pressure correction a mixture pressure 127 ! correction equation by adding solids phases that do not become 128 ! close packed 129 DO M = 1, MMAX 130 IF (.NOT.CLOSE_PACKED(M)) THEN 131 DO IJK = ijkstart3, ijkend3 132 IF (FLUID_AT(IJK)) THEN 133 IPJK = IP_OF(IJK) 134 IJPK = JP_OF(IJK) 135 IJKP = KP_OF(IJK) 136 137 ! East face (i+1/2, j, k) 138 AM = ROP_SE(IJK,M)*AYZ(IJK) 139 A_M(IJK,E,M) = AM 140 A_M(IPJK,W,M) = AM 141 142 ! North face (i, j+1/2, k) 143 AM = ROP_SN(IJK,M)*AXZ(IJK) 144 A_M(IJK,N,M) = AM 145 A_M(IJPK,S,M) = AM 146 147 ! Top face (i, j, k+1/2) 148 IF (DO_K) THEN 149 AM = ROP_ST(IJK,M)*AXY(IJK) 150 A_M(IJK,T,M) = AM 151 A_M(IJKP,B,M) = AM 152 ENDIF 153 154 ! West face (i-1/2, j, k) 155 IMJK = IM_OF(IJK) 156 IF (.NOT.FLUID_AT(IMJK)) THEN 157 AM = ROP_SE(IMJK,M)*AYZ(IMJK) 158 A_M(IJK,W,M) = AM 159 ENDIF 160 161 ! South face (i, j-1/2, k) 162 IJMK = JM_OF(IJK) 163 IF (.NOT.FLUID_AT(IJMK)) THEN 164 AM = ROP_SN(IJMK,M)*AXZ(IJMK) 165 A_M(IJK,S,M) = AM 166 ENDIF 167 168 ! Bottom face (i, j, k-1/2) 169 IF (DO_K) THEN 170 IJKM = KM_OF(IJK) 171 IF (.NOT.FLUID_AT(IJKM)) THEN 172 AM = ROP_ST(IJKM,M)*AXY(IJKM) 173 A_M(IJK,B,M) = AM 174 ENDIF 175 ENDIF 176 177 ENDIF ! end if(fluid_at(ijk)) 178 ENDDO ! end do ijk=ijkstart3,ijkend3 179 ENDIF ! end if(.not.close_packed) 180 ENDDO ! end do (m=1,mmax) 181 182 RETURN 183 END SUBROUTINE CONV_PP_G 184