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 calcgammak(nface,ielfa,vel,gradkel,gamma,xlet,
20     &  xxn,xxj,ipnei,betam,nef,flux)
21!
22!     determine gamma for the turbulent kinetic energy:
23!        upwind difference: gamma=0
24!        central difference: gamma=1
25!
26      implicit none
27!
28      integer nface,ielfa(4,*),i,j,indexf,ipnei(*),iel1,iel2,nef
29!
30      real*8 vel(nef,0:7),gradkel(3,*),xxn(3,*),xxj(3,*),vud,vcd,
31     &  gamma(*),phic,xlet(*),betam,flux(*)
32!
33      do i=1,nface
34         iel2=ielfa(2,i)
35!
36!        faces with only one neighbor need not be treated
37!
38         if(iel2.le.0) cycle
39         iel1=ielfa(1,i)
40         j=ielfa(4,i)
41         indexf=ipnei(iel1)+j
42!
43         vcd=vel(iel2,6)-vel(iel1,6)
44         if(dabs(vcd).lt.1.d-3*dabs(vel(iel1,6))) vcd=0.d0
45!
46         if(flux(indexf).ge.0.d0) then
47            vud=2.d0*xlet(indexf)*
48     &           (gradkel(1,iel1)*xxj(1,indexf)+
49     &            gradkel(2,iel1)*xxj(2,indexf)+
50     &            gradkel(3,iel1)*xxj(3,indexf))
51         else
52            vud=2.d0*xlet(indexf)*
53     &           (gradkel(1,iel2)*xxj(1,indexf)+
54     &            gradkel(2,iel2)*xxj(2,indexf)+
55     &            gradkel(3,iel2)*xxj(3,indexf))
56         endif
57!
58         if(dabs(vud).lt.1.d-20) then
59            gamma(i)=0.d0
60            cycle
61         endif
62!
63         phic=1.d0-vcd/vud
64!
65         if(phic.ge.1.d0) then
66            gamma(i)=0.d0
67         elseif(phic.le.0.d0) then
68            gamma(i)=0.d0
69         elseif(betam.le.phic) then
70            gamma(i)=1.d0
71         else
72            gamma(i)=phic/betam
73         endif
74c         write(*,*) 'calcgammat',i,gamma(i)
75      enddo
76!
77      return
78      end
79