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!     This program is distributed in the hope that it will be useful,
10!     but WITHOUT ANY WARRANTY; without even the implied warranty of
11!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!     GNU General Public License for more details.
13!
14!     You should have received a copy of the GNU General Public License
15!     along with this program; if not, write to the Free Software
16!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17!
18!     author: Yannick Muller
19!
20      subroutine calc_ider_wye(df,pt1,Tt1,xflow1,xflow2,pt2,
21     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag
22     &,zeta)
23!
24      implicit none
25!
26      integer ichan_num,ider,iflag
27!
28      real*8
29     &df(6),
30     &pt1,
31     &pt2,
32     &Tt1,
33     &Tt2,
34     &xflow1,
35     &xflow2,
36     &A1,
37     &A2,
38     &A_s,
39     &kappa,
40     &R,
41     &dh1,
42     &dh2,
43     &alpha,
44     &calc_residual_wye,
45     &eps,
46     &h,
47     &f0,
48     &zeta_fac,
49     &zeta
50!
51      eps = 1.0e-4
52!
53      f0 = calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2,
54     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,
55     &R,ider,iflag,zeta)
56!
57      h = eps*dabs(pt1)
58      if(h.eq.0)then
59         h = eps
60      endif
61      df(1) = (calc_residual_wye(pt1+h,Tt1,xflow1,xflow2,pt2,
62     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
63     &iflag,zeta)-f0)/h
64!
65      h = eps*dabs(Tt1)
66      if(h.eq.0)then
67         h = eps
68      endif
69      df(2) = (calc_residual_wye(pt1,Tt1+h,xflow1,xflow2,pt2,
70     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
71     &iflag,zeta)-f0)/h
72!
73      h = eps*dabs(xflow1)
74      if(h.eq.0)then
75         h = eps
76      endif
77      df(3) = (calc_residual_wye(pt1,Tt1,xflow1+h,xflow2,pt2,
78     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
79     &iflag,zeta)-f0)/h
80!
81      h = eps*dabs(xflow2)
82      if(h.eq.0)then
83         h = eps
84      endif
85      df(4) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2+h,pt2,
86     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
87     &iflag,zeta)-f0)/h
88!
89      h = eps*dabs(pt2)
90      if(h.eq.0)then
91         h = eps
92      endif
93      df(5) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2+h,
94     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
95     &iflag,zeta)-f0)/h
96!
97      h = eps*dabs(Tt2)
98      if(h.eq.0)then
99         h = eps
100      endif
101      df(6) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2,
102     &Tt2+h,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,
103     &iflag,zeta)-f0)/h
104!
105      return
106      end
107