1 subroutine genBC1 (BCtmp, iBC, BC) 2c 3c---------------------------------------------------------------------- 4c 5c This subroutine adjusts the boundary conditions to accommodate for 6c the velocity constraints in the non-axes directions. It copies the 7c reduced constraint parameters in BC. 8c 9c input: 10c BCtmp (nshg,6+5*I3nsd) : input BC parameters (density, temperature, 11c pressure, (nsd-1)(nsd+1) velocity params, 12c upto 4 scalar params) 13c iBC (nshg) : boundary condition code 14c 15c output: 16c BC (nshg,ndofBC) : the constraint eq's parameters 17c 18c 19c Farzin Shakib, Winter 1986. 20c Zdenek Johan, Winter 1991. (Fortran 90) 21c---------------------------------------------------------------------- 22c 23 include "common.h" 24 25c 26 dimension BCtmp(nshg,ndof+7), iBC(nshg), 27 & BC(nshg,ndofBC),tmpbc(4) 28c 29 dimension tmp(nshg) 30c 31c.... scalars 32c 33 do isclr=1,nsclr 34 where (btest(iBC,5+isclr)) BC(:,6+isclr) = BCtmp(:,12+isclr) 35 enddo 36c 37c.... set up the thermodynamic properties 38c 39 where (btest(iBC,0)) BC(:,1) = BCtmp(:,1) ! density 40 where (btest(iBC,1)) BC(:,2) = BCtmp(:,2) ! temperature 41 where (btest(iBC,2)) BC(:,1) = BCtmp(:,3) ! pressure 42c 43c.... if the velocity in the x1-direction is specified 44c 45 where (ibits(iBC,3,3) .eq. 1) 46 tmp = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2 47 BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,4) 48 BC(:,4) = BCtmp(:,5) / BCtmp(:,4) 49 BC(:,5) = BCtmp(:,6) / BCtmp(:,4) 50 endwhere 51c 52c.... if the velocity in the x2-direction is specified 53c 54 where (ibits(iBC,3,3) .eq. 2) 55 tmp = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2 56 BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,5) 57 BC(:,4) = BCtmp(:,4) / BCtmp(:,5) 58 BC(:,5) = BCtmp(:,6) / BCtmp(:,5) 59 endwhere 60c 61c.... if the two velocities are specified (x1 & x2-direction) 62c 63c 64c Protect against user flipping the order of x1 and x2 in 65c the vector 1 and vector 2. Without this it will blow up. 66c 67 do i=1,nshg 68 if(ibits(iBC(i),3,3) .eq. 3 .and. 69 & (BCtmp(i,4).eq.0 .or. BCtmp(i,9).eq.0)) then !flip them 70 tmpbc(1:4)=BCtmp(i,4:7) 71 BCtmp(i,4:7)=BCtmp(i,8:11) 72 BCtmp(i,8:11)=tmpbc(1:4) 73 endif 74 enddo 75 where (ibits(iBC,3,3) .eq. 3) 76 tmp = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2 77 & + BCtmp(:, 6)**2) 78 BCtmp(:, 4) = BCtmp(:, 4) / tmp 79 BCtmp(:, 5) = BCtmp(:, 5) / tmp 80 BCtmp(:, 6) = BCtmp(:, 6) / tmp 81 BCtmp(:, 7) = BCtmp(:, 7) * tmp 82c 83 tmp = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2 84 & + BCtmp(:,10)**2) 85 BCtmp(:, 8) = BCtmp(:, 8) / tmp 86 BCtmp(:, 9) = BCtmp(:, 9) / tmp 87 BCtmp(:,10) = BCtmp(:,10) / tmp 88 BCtmp(:,11) = BCtmp(:,11) * tmp 89c 90 BCtmp(:, 4) = BCtmp(:, 9) * BCtmp(:, 4) 91 & - BCtmp(:, 5) * BCtmp(:, 8) 92 BCtmp(:, 6) = BCtmp(:, 9) * BCtmp(:, 6) 93 & - BCtmp(:, 5) * BCtmp(:,10) 94 BCtmp(:, 7) = BCtmp(:, 9) * BCtmp(:, 7) 95 & - BCtmp(:, 5) * BCtmp(:,11) 96 BC(:,3) = BCtmp(:, 7) / BCtmp(:, 4) 97 BC(:,4) = BCtmp(:, 6) / BCtmp(:, 4) 98c 99 BCtmp(:, 9) = BCtmp(:, 4) * BCtmp(:, 9) 100 BCtmp(:,10) = BCtmp(:, 4) * BCtmp(:,10) 101 & - BCtmp(:, 8) * BCtmp(:, 6) 102 BCtmp(:,11) = BCtmp(:, 4) * BCtmp(:,11) 103 & - BCtmp(:, 8) * BCtmp(:, 7) 104 BC(:,5) = BCtmp(:,11) / BCtmp(:, 9) 105 BC(:,6) = BCtmp(:,10) / BCtmp(:, 9) 106 endwhere 107c 108c.... if the velocity in the x3-direction is specified 109c 110 if (nsd .eq. 3) then 111 where (ibits(iBC,3,3) .eq. 4) 112 tmp = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2 113 BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,6) 114 BC(:,4) = BCtmp(:,4) / BCtmp(:,6) 115 BC(:,5) = BCtmp(:,5) / BCtmp(:,6) 116 endwhere 117 endif 118c 119c.... if two velocities are specified (x1 & x3-direction) 120c 121 if (nsd .eq. 3) then 122c 123c Protect against user flipping the order of x1 and x3 in 124c the vector 1 and vector 2. Without this it will blow up. 125c 126 do i=1,nshg 127 if(ibits(iBC(i),3,3) .eq. 5 .and. 128 & (BCtmp(i,4).eq.0 .or. BCtmp(i,10).eq.0)) then !flip them 129 tmpbc(1:4)=BCtmp(i,4:7) 130 BCtmp(i,4:7)=BCtmp(i,8:11) 131 BCtmp(i,8:11)=tmpbc(1:4) 132 endif 133 enddo 134 where (ibits(iBC,3,3) .eq. 5) 135 tmp = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2 136 & + BCtmp(:, 6)**2) 137 BCtmp(:, 4) = BCtmp(:, 4) / tmp 138 BCtmp(:, 5) = BCtmp(:, 5) / tmp 139 BCtmp(:, 6) = BCtmp(:, 6) / tmp 140 BCtmp(:, 7) = BCtmp(:, 7) * tmp 141c 142 tmp = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2 143 & + BCtmp(:,10)**2) 144 BCtmp(:, 8) = BCtmp(:, 8) / tmp 145 BCtmp(:, 9) = BCtmp(:, 9) / tmp 146 BCtmp(:,10) = BCtmp(:,10) / tmp 147 BCtmp(:,11) = BCtmp(:,11) * tmp 148c 149 BCtmp(:, 4) = BCtmp(:,10) * BCtmp(:, 4) 150 & - BCtmp(:, 6) * BCtmp(:, 8) 151 BCtmp(:, 5) = BCtmp(:,10) * BCtmp(:, 5) 152 & - BCtmp(:, 6) * BCtmp(:, 9) 153 BCtmp(:, 7) = BCtmp(:,10) * BCtmp(:, 7) 154 & - BCtmp(:, 6) * BCtmp(:,11) 155 BC(:,3) = BCtmp(:, 7) / BCtmp(:, 4) 156 BC(:,4) = BCtmp(:, 5) / BCtmp(:, 4) 157c 158 BCtmp(:, 9) = BCtmp(:, 4) * BCtmp(:, 9) 159 & - BCtmp(:, 8) * BCtmp(:, 5) 160 BCtmp(:,10) = BCtmp(:, 4) * BCtmp(:,10) 161 BCtmp(:,11) = BCtmp(:, 4) * BCtmp(:,11) 162 & - BCtmp(:, 8) * BCtmp(:, 7) 163 BC(:,5) = BCtmp(:,11) / BCtmp(:,10) 164 BC(:,6) = BCtmp(:, 9) / BCtmp(:,10) 165 endwhere 166 endif 167c 168c.... if two velocities are specified (x2 & x3-direction) 169c 170 if (nsd .eq. 3) then 171c 172c Protect against user flipping the order of x2 and x3 in 173c the vector 1 and vector 2. Without this it will blow up. 174c 175 do i=1,nshg 176 if(ibits(iBC(i),3,3) .eq. 6 .and. ( 177 & BCtmp(i,5).eq.0 .or. BCtmp(i,10).eq.0)) then !flip them 178 tmpbc(1:4)=BCtmp(i,4:7) 179 BCtmp(i,4:7)=BCtmp(i,8:11) 180 BCtmp(i,8:11)=tmpbc(1:4) 181 endif 182 enddo 183 where (ibits(iBC,3,3) .eq. 6) 184 tmp = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2 185 & + BCtmp(:, 6)**2) 186 BCtmp(:, 4) = BCtmp(:, 4) / tmp 187 BCtmp(:, 5) = BCtmp(:, 5) / tmp 188 BCtmp(:, 6) = BCtmp(:, 6) / tmp 189 BCtmp(:, 7) = BCtmp(:, 7) * tmp 190c 191 tmp = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2 192 & + BCtmp(:,10)**2) 193 BCtmp(:, 8) = BCtmp(:, 8) / tmp 194 BCtmp(:, 9) = BCtmp(:, 9) / tmp 195 BCtmp(:,10) = BCtmp(:,10) / tmp 196 BCtmp(:,11) = BCtmp(:,11) * tmp 197c 198 BCtmp(:, 4) = BCtmp(:,10) * BCtmp(:, 4) 199 & - BCtmp(:, 6) * BCtmp(:, 8) 200 BCtmp(:, 5) = BCtmp(:,10) * BCtmp(:, 5) 201 & - BCtmp(:, 6) * BCtmp(:, 9) 202 BCtmp(:, 7) = BCtmp(:,10) * BCtmp(:, 7) 203 & - BCtmp(:, 6) * BCtmp(:,11) 204 BC(:,3) = BCtmp(:, 7) / BCtmp(:, 5) 205 BC(:,4) = BCtmp(:, 4) / BCtmp(:, 5) 206c 207 BCtmp(:, 8) = BCtmp(:, 5) * BCtmp(:, 8) 208 & - BCtmp(:, 9) * BCtmp(:, 4) 209 BCtmp(:,10) = BCtmp(:, 5) * BCtmp(:,10) 210 BCtmp(:,11) = BCtmp(:, 5) * BCtmp(:,11) 211 & - BCtmp(:, 9) * BCtmp(:, 7) 212 BC(:,5) = BCtmp(:,11) / BCtmp(:,10) 213 BC(:,6) = BCtmp(:, 8) / BCtmp(:,10) 214 endwhere 215 endif 216c 217c.... if all velocities are specified 218c 219 if (nsd .eq. 3) then 220 where (ibits(iBC,3,3) .eq. 7) 221 BC(:,3) = BCtmp(:,7) * BCtmp(:,4) 222 BC(:,4) = BCtmp(:,7) * BCtmp(:,5) 223 BC(:,5) = BCtmp(:,7) * BCtmp(:,6) 224 endwhere 225 endif 226c 227c.... end 228c 229 return 230 end 231