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_tee(df,pt1,Tt1,xflow1,xflow2,pt2,
21     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)
22!
23      implicit none
24!
25      integer ider,iflag
26!
27      real*8
28     &df(6),
29     &pt1,
30     &pt2,
31     &Tt1,
32     &Tt2,
33     &xflow1,
34     &xflow2,
35     &A1,
36     &A2,
37     &kappa,
38     &R,
39     &calc_residual_tee,
40!     Accuracy of the numerical derivatives
41     &eps,
42!     Step for the derivative
43     &h,
44     &f0,
45     &zeta_fac,
46     &zeta
47!
48      eps = 1.0e-8
49!
50      f0 = calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2,
51     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)
52!
53      h = eps*dabs(pt1)
54      if(h.eq.0)then
55         h = eps
56      endif
57      df(1) = (calc_residual_tee(pt1+h,Tt1,xflow1,xflow2,pt2,
58     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
59!
60      h = eps*dabs(Tt1)
61      if(h.eq.0)then
62         h = eps
63      endif
64      df(2) = (calc_residual_tee(pt1,Tt1+h,xflow1,xflow2,pt2,
65     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
66!
67      h = eps*dabs(xflow1)
68      if(h.eq.0)then
69         h = eps
70      endif
71      df(3) = (calc_residual_tee(pt1,Tt1,xflow1+h,xflow2,pt2,
72     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
73!
74      h = eps*dabs(xflow2)
75      if(h.eq.0)then
76         h = eps
77      endif
78      df(4) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2+h,pt2,
79     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
80!
81      h = eps*dabs(pt2)
82      if(h.eq.0)then
83         h = eps
84      endif
85      df(5) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2+h,
86     &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
87!
88      h = eps*dabs(Tt2)
89      if(h.eq.0)then
90         h = eps
91      endif
92      df(6) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2,
93     &Tt2+h,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h
94!
95      return
96      end
97