1 subroutine e3qvar (ycl, shp, shgl, rho, 2 & xl, g1yi, g2yi, g3yi, 3 & shg, dxidx, WdetJ, T, 4 & cp, u1, u2, u3) 5c 6c---------------------------------------------------------------------- 7c 8c This routine computes the variables at integration point 9c necessary for the computation of the diffusive flux vector. 10c 11c input: 12c ycl (npro,nshape,ndof) : primitive variables 13c shp (npro,nshape) : element shape-functions 14c shgl (npro,nsd,nshape) : element local-grad-shape-functions 15c xl (npro,nenl,nsd) : nodal coordinates at current step 16c 17c output: 18c g1yi (npro,nflow) : grad-y in direction 1 19c g2yi (npro,nflow) : grad-y in direction 2 20c g3yi (npro,nflow) : grad-y in direction 3 21c shg (npro,nshape,nsd) : element global grad-shape-functions 22c dxidx (npro,nsd,nsd) : inverse of deformation gradient 23c WdetJ (npro) : weighted Jacobian 24c T (npro) : temperature 25c cp (npro) : specific heat at constant pressure 26c u1 (npro) : x1-velocity component 27c u2 (npro) : x2-velocity component 28c u3 (npro) : x3-velocity component 29c 30c 31c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 32c Zdenek Johan, Winter 1991. (Fortran 90) 33c Kenneth Jansen, Winter 1997. Primitive Variables 34c---------------------------------------------------------------------- 35c 36 include "common.h" 37c 38c passed arrays 39c 40 dimension ycl(npro,nshl,ndof), rho(npro), 41 & shp(npro,nshl), tmp(npro), 42 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 43 & g1yi(npro,nflow), g2yi(npro,nflow), 44 & g3yi(npro,nflow), shg(npro,nshl,nsd), 45 & dxidx(npro,nsd,nsd), WdetJ(npro), 46 & T(npro), cp(npro), 47 & u1(npro), u2(npro), 48 & u3(npro), pres(npro) 49c 50c local arrays 51c 52 dimension dxdxi(npro,nsd,nsd), Sclr(npro) 53 54 ttim(34) = ttim(34) - secs(0.0) 55c 56c.... -------------> Primitive variables at int. point <-------------- 57c 58c.... compute primitive variables 59c 60 pres = zero 61 u1 = zero 62 u2 = zero 63 u3 = zero 64 T = zero 65 66 67c 68 do n = 1, nshl 69c 70c y(intp)=SUM_{a=1}^nshape (N_a(intp) Ya) (we don't need pressure here) 71c 72 pres = pres + shp(:,n) * ycl(:,n,1) 73 u1 = u1 + shp(:,n) * ycl(:,n,2) 74 u2 = u2 + shp(:,n) * ycl(:,n,3) 75 u3 = u3 + shp(:,n) * ycl(:,n,4) 76 T = T + shp(:,n) * ycl(:,n,5) 77 enddo 78c 79c$$$c.... compute cp, only thermodynamic property needed for qdiff 80c$$$c 81c$$$ cp = Rgas * gamma / gamma1 82c.... get the thermodynamic properties 83c 84 85 if (iLSet .ne. 0)then 86 Sclr = zero 87 isc=abs(iRANS)+6 88 do n = 1, nshl 89 Sclr = Sclr + shp(:,n) * ycl(:,n,isc) 90 enddo 91 endif 92c 93 if (Navier .eq. 1) ithm = 7 94 95 call getthm (pres, T, Sclr, 96 & tmp, rho, tmp, 97 & tmp, tmp, tmp, 98 & cp, tmp, tmp, 99 & tmp, tmp) 100c 101c.... ---------------------> Element Metrics <----------------------- 102c 103 call e3metric( xl, shgl, dxidx, 104 & shg, WdetJ) 105 106c 107c.... ---------------------> Global Gradients <----------------------- 108c 109 g1yi = zero 110 g2yi = zero 111 g3yi = zero 112c 113c 114 do n = 1, nshl 115c.... compute the global gradient of Y-variables 116c 117c 118c Y_{,x_i}=SUM_{a=1}^nenl (N_{a,x_i}(intp) Ya) 119c 120 g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * ycl(:,n,1) 121 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * ycl(:,n,2) 122 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * ycl(:,n,3) 123 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * ycl(:,n,4) 124 g1yi(:,5) = g1yi(:,5) + shg(:,n,1) * ycl(:,n,5) 125c 126 g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * ycl(:,n,1) 127 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * ycl(:,n,2) 128 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * ycl(:,n,3) 129 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * ycl(:,n,4) 130 g2yi(:,5) = g2yi(:,5) + shg(:,n,2) * ycl(:,n,5) 131c 132 g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * ycl(:,n,1) 133 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * ycl(:,n,2) 134 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * ycl(:,n,3) 135 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * ycl(:,n,4) 136 g3yi(:,5) = g3yi(:,5) + shg(:,n,3) * ycl(:,n,5) 137 138 enddo 139 140 ttim(34) = ttim(34) + secs(0.0) 141 142c 143c.... return 144c 145 146 return 147 end 148