1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine calcgammav(ielfa,vel,gradvel,gamma,xlet,
20     &  xxj,ipnei,betam,nef,flux,nfacea,nfaceb)
21!
22!     determine gamma for the velocity:
23!        upwind difference: gamma=0
24!        central difference: gamma=1
25!
26      implicit none
27!
28      integer ielfa(4,*),i,j,indexf,ipnei(*),iel1,iel2,nef,
29     &  nfacea,nfaceb
30!
31      real*8 vel(nef,0:7),gradvel(3,3,*),xxj(3,*),vud,vcd,
32     &  gamma(*),phic,xlet(*),betam,flux(*),dvel1,dvel2
33!
34!
35!
36      do i=nfacea,nfaceb
37         iel2=ielfa(2,i)
38!
39!        faces with only one neighbor need not be treated
40!
41         if(iel2.le.0) cycle
42         iel1=ielfa(1,i)
43         j=ielfa(4,i)
44         indexf=ipnei(iel1)+j
45!
46         dvel1=dsqrt(vel(iel1,1)**2+vel(iel1,2)**2+vel(iel1,3)**2)
47         dvel2=dsqrt(vel(iel2,1)**2+vel(iel2,2)**2+vel(iel2,3)**2)
48!
49         vcd=dvel2-dvel1
50!
51         if(dabs(vcd).lt.1.d-3*dvel1) vcd=0.d0
52!
53         if(flux(indexf).ge.0.d0) then
54!
55            vud=2.d0*xlet(indexf)*(
56     &       (vel(iel1,1)*gradvel(1,1,iel1)+
57     &        vel(iel1,2)*gradvel(2,1,iel1)+
58     &        vel(iel1,3)*gradvel(3,1,iel1))*xxj(1,indexf)+
59     &       (vel(iel1,1)*gradvel(1,2,iel1)+
60     &        vel(iel1,2)*gradvel(2,2,iel1)+
61     &        vel(iel1,3)*gradvel(3,2,iel1))*xxj(2,indexf)+
62     &       (vel(iel1,1)*gradvel(1,3,iel1)+
63     &        vel(iel1,2)*gradvel(2,3,iel1)+
64     &        vel(iel1,3)*gradvel(3,3,iel1))*xxj(3,indexf))
65            vcd=vcd*dvel1
66         else
67            vud=2.d0*xlet(indexf)*(
68     &       (vel(iel2,1)*gradvel(1,1,iel2)+
69     &        vel(iel2,2)*gradvel(2,1,iel2)+
70     &        vel(iel2,3)*gradvel(3,1,iel2))*xxj(1,indexf)+
71     &       (vel(iel2,1)*gradvel(1,2,iel2)+
72     &        vel(iel2,2)*gradvel(2,2,iel2)+
73     &        vel(iel2,3)*gradvel(3,2,iel2))*xxj(2,indexf)+
74     &       (vel(iel2,1)*gradvel(1,3,iel2)+
75     &        vel(iel2,2)*gradvel(2,3,iel2)+
76     &        vel(iel2,3)*gradvel(3,3,iel2))*xxj(3,indexf))
77            vcd=vcd*dvel2
78         endif
79!
80         if(dabs(vud).lt.1.d-20) then
81            gamma(i)=0.d0
82            cycle
83         endif
84!
85         phic=1.d0-vcd/vud
86!
87         if(phic.ge.1.d0) then
88            gamma(i)=0.d0
89         elseif(phic.le.0.d0) then
90            gamma(i)=0.d0
91         elseif(betam.le.phic) then
92            gamma(i)=1.d0
93         else
94            gamma(i)=phic/betam
95         endif
96      enddo
97!
98      return
99      end
100