1 subroutine itrBC (y,ac, iBC, BC, iper, ilwork) 2c 3c---------------------------------------------------------------------- 4c 5c This program satisfies the boundary conditions on the Y-variables. 6c 7c input: 8c y (nshg,nflow) : y variables 9c iBC (nshg) : Boundary Condition Code 10c BC (nshg,ndofBC) : boundary condition constraint parameters 11c ylimit (3,nflow) : (1,:) limiting flag 12c (2,:) lower bound 13c (3,:) upper bound 14c output: 15c y (nshg,nflow) : Adjusted V value(s) corresponding to a 16c constraint d.o.f. 17c 18c 19c Farzin Shakib, Winter 1987. 20c Zdenek Johan, Winter 1991. (Fortran 90) 21c---------------------------------------------------------------------- 22c 23 include "common.h" 24c 25 dimension y(nshg,nflow), iBC(nshg), 26 & ac(nshg,nflow), BC(nshg,ndofBC) 27 28 dimension ilwork(nlwork), iper(nshg) 29 30 real*8 tmp(nshg), y1(nshg),q1(nshg) 31 dimension rn1(nshg), rmagn1(nshg) 32 real*8 limitcount(nflow) 33 integer locmax(1),locmin(1) 34c 35c limiting...ugly but sometimes the only way 36c 37 limitcount=0 38 do i=1,nflow 39 if(ylimit(1,i).gt.0) then 40 locmax=maxloc(y(:,i)) 41 locmin=minloc(y(:,i)) 42 ymax=maxval(y(:,i)) 43 ymin=minval(y(:,i)) 44c write(33,34)i,ymax,ymin,1.*locmax,1.*locmin 45c call flush(33) 46 do j=1,numnp ! only limit linear modes for now 47 ypc=y(j,i) 48 y(j,i)=min(ylimit(3,i),max(ylimit(2,i),y(j,i))) 49 if(ypc.ne.y(j,i) )limitcount(i)=limitcount(i)+one 50 enddo 51 endif 52 enddo 53c if(maxval(limitcount).gt.0) then 54c write(33,34)myrank,(limitcount(j)/numnp,j=1,nflow) 55c call flush(33) 56c endif 57 34 format(i5,5(2x,e14.7)) 58c 59c.... ------------------------> Temperature <------------------------ 60c.... temperature 61c 62 where (btest(iBC,1)) 63 y(:,5) = BC(:,2) 64 endwhere 65c 66c.... -------------------------> Velocity <-------------------------- 67c.... 3D 68c 69c.... x1-velocity, 3D 70c 71 where (ibits(iBC,3,3) .eq. 1) 72c$$$ rn1(:)=sqrt(one/(one + BC(:,4)**2+BC(:,4)**2)) 73c$$$ rmagn1(:)=rn1(:)**2*(y(:,1)+y(:,2)*BC(:,4) 74c$$$ & +y(:,3)*BC(:,5)-BC(:,3)) 75c$$$ y(:,1) = y(:,1)-rmagn1(:) 76c$$$ y(:,2) = y(:,2)-rmagn1(:)*BC(:,4) 77c$$$ y(:,3) = y(:,3)-rmagn1(:)*BC(:,5) 78 79 y(:,1) = BC(:,3) - BC(:,4) * y(:,2) 80 & - BC(:,5) * y(:,3) 81 endwhere 82c 83c.... x2-velocity, 3D 84c 85 where (ibits(iBC,3,3) .eq. 2) 86 y(:,2) = BC(:,3) - BC(:,4) * y(:,1) 87 & - BC(:,5) * y(:,3) 88 endwhere 89c 90c.... x1-velocity and x2-velocity, 3D 91c 92 where (ibits(iBC,3,3) .eq. 3) 93 y(:,1) = BC(:,3) - BC(:,4) * y(:,3) 94 y(:,2) = BC(:,5) - BC(:,6) * y(:,3) 95 endwhere 96c 97c.... x3-velocity, 3D 98c 99 where (ibits(iBC,3,3) .eq. 4) 100 y(:,3) = BC(:,3) - BC(:,4) * y(:,1) 101 & - BC(:,5) * y(:,2) 102 endwhere 103c 104c.... x1-velocity and x3-velocity, 3D 105c 106 where (ibits(iBC,3,3) .eq. 5) 107 y(:,1) = BC(:,3) - BC(:,4) * y(:,2) 108 y(:,3) = BC(:,5) - BC(:,6) * y(:,2) 109 endwhere 110c 111c.... x2-velocity and x3-velocity, 3D 112c 113 where (ibits(iBC,3,3) .eq. 6) 114 y(:,2) = BC(:,3) - BC(:,4) * y(:,1) 115 y(:,3) = BC(:,5) - BC(:,6) * y(:,1) 116 endwhere 117c 118c.... x1-velocity, x2-velocity and x3-velocity, 3D 119c 120 where (ibits(iBC,3,3) .eq. 7) 121 y(:,1) = BC(:,3) 122 y(:,2) = BC(:,4) 123 y(:,3) = BC(:,5) 124 endwhere 125c 126c endif 127c 128c.... end of velocity 129c 130c.... --------------------------> Density <-------------------------- 131c 132 if (any(btest(iBC,0))) then 133c 134c.... density 135c 136 where (btest(iBC,0)) 137 q1 = BC(:,1) ! density 138 elsewhere 139 q1 = one 140 endwhere 141c 142c 143c 144 npro = nshg 145c 146 ithm = 2 ! get pressure from rho and T 147c...when ithm=2 scalar is not used so tmp is in place 148 call getthm (y1, y(:,5), tmp, 149 & rk, q1, tmp, 150 & tmp, tmp, tmp, 151 & tmp, tmp, tmp, 152 & tmp, tmp) 153c 154 where (btest(iBC,0)) 155 y(:,1) = y1 156 endwhere 157 158c 159 endif 160c 161c.... -------------------------> Pressure <-------------------------- 162c 163 if (any(btest(iBC,2))) then 164c 165c.... pressure 166c 167 where (btest(iBC,2)) 168 y(:,4) = BC(:,1) ! pressure here 169 endwhere 170c 171 endif 172c 173c.... local periodic (and axisymmetric) boundary conditions (no communications) 174c 175 do i = 1,nflow 176 y(:,i) = y(iper(:),i) 177 if(ires.ne.2) ac(:,i) = ac(iper(:),i) 178 enddo 179c 180c.... communications 181c 182 if (numpe > 1) then 183 call commu (y, ilwork, nflow, 'out') 184 if(ires.ne.2) call commu (ac, ilwork, nflow, 'out') 185 endif 186c 187c slave has masters value, for abc we need to rotate it 188c 189 if(iabc==1) then !are there any axisym bc's 190 call rotabc(y, iBC, 'out') 191 if(ires.ne.2) call rotabc(ac, iBC, 'out') 192 endif 193 194c 195c.... return 196c 197 return 198 end 199 200 201 subroutine itrBCSclr (y, ac, iBC, BC, iper, ilwork) 202c 203c---------------------------------------------------------------------- 204c 205c This routine satisfies the boundary conditions on the isclr 206c 207c---------------------------------------------------------------------- 208c 209 include "common.h" 210c 211 dimension y(nshg,ndof), iBC(nshg), 212 & ac(nshg,ndof), BC(nshg,ndofBC) 213 214 dimension ilwork(nlwork), iper(nshg) 215 dimension T(nshg) 216 217 if(isclr.ne.0) then 218 id=5+isclr 219 ibb=6+isclr 220 ib=4+isclr 221 endif 222c 223c limiting...ugly but sometimes the only way 224c 225 if(ylimit(1,id).gt.0) 226 & y(:,id)=min(ylimit(3,id),max(ylimit(2,id),y(:,id))) 227c 228c.... ------------------------> Scalar <------------------------ 229c 230c 231 where (btest(iBC,id)) 232 y(:,id) = BC(:,ibb) 233 endwhere 234c 235c.... local periodic (and axisymmetric) boundary conditions (no communications) 236c 237 do i = 1,nshg 238c if (btest(iBC(i),10)) then 239 y(i,id) = y(iper(i),id) 240 ac(i,id) = ac(iper(i),id) 241c end if 242 enddo 243c 244c.... communications 245c 246 if (numpe > 1) then 247 T=y(:,id) 248 call commu (T, ilwork, 1, 'out') 249 y(:,id)=T 250 T=ac(:,id) 251 call commu (T, ilwork, 1, 'out') 252 ac(:,id)=T 253 endif 254 255 ttim(53) = ttim(53) + tmr() 256c 257 return 258 end 259 260 261 262