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 correctflux(nef,ipnei,neifa,neiel,flux,vfa,advfa,area,
20     &  vel,alet,ielfa,ale,ifabou,nefa,nefb,xxnj,gradpcfa)
21!
22!     correction of v due to the balance of mass
23!     the correction is in normal direction to the face
24!
25      implicit none
26!
27      integer i,nef,indexf,ipnei(*),neifa(*),neiel(*),ielfa(4,*),
28     &  iel,ifa,ifabou(*),indexb,nefa,nefb
29!
30      real*8 flux(*),vfa(0:7,*),advfa(*),area(*),vel(nef,0:7),alet(*),
31     &  ale(*),xxnj(3,*),gradpcfa(3,*)
32!
33!
34!
35      do i=nefa,nefb
36c         totflux=0.d0
37         do indexf=ipnei(i)+1,ipnei(i+1)
38            ifa=neifa(indexf)
39            iel=neiel(indexf)
40            if(iel.gt.0) then
41!
42!              internal face
43!
44               flux(indexf)=flux(indexf)+vfa(5,ifa)*advfa(ifa)*
45     &                      ((vel(i,4)-vel(iel,4))*alet(indexf)
46     &                 -(gradpcfa(1,ifa)*xxnj(1,indexf)+
47     &                   gradpcfa(2,ifa)*xxnj(2,indexf)+
48     &                   gradpcfa(3,ifa)*xxnj(3,indexf)))
49            else
50               indexb=-ielfa(2,ifa)
51               if(indexb.gt.0) then
52                  if(((ifabou(indexb+1).eq.0).or.
53     &                 (ifabou(indexb+2).eq.0).or.
54     &                 (ifabou(indexb+3).eq.0)).and.
55     &                 (ifabou(indexb+4).ne.0)) then
56!
57!                    external face with pressure boundary conditions
58!
59                     flux(indexf)=flux(indexf)
60     &                      +vfa(5,ifa)*advfa(ifa)*(vel(i,4)*ale(indexf)
61     &                             -(gradpcfa(1,ifa)*xxnj(1,indexf)+
62     &                               gradpcfa(2,ifa)*xxnj(2,indexf)+
63     &                               gradpcfa(3,ifa)*xxnj(3,indexf)))
64                  endif
65               endif
66            endif
67c            totflux=totflux+flux(indexf)
68         enddo
69      enddo
70!
71      return
72      end
73