1!     CalculiX - A 3-dimensional finite element program
2!              Copyright (C) 1998-2021 Guido Dhondt
3!
4!     This program is free software; you can redistribute it and/or
5!     modify it under the terms of the GNU General Public License as
6!     published by the Free Software Foundation(version 2);
7!
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      subroutine materialdata_cp_sec(imat,ntmat_,t1l,shcon,nshcon,cp,
19     &  physcon)
20!
21      implicit none
22!
23!     determines the secant specific heat at constant pressure cp
24!
25!     the difference with materialdata_cp is that the specific heat at
26!     constant pressure cp as returned from the present routine
27!     is the secant value and not the differential value.
28!     For the differential value we have:
29!            dh=cp*dT
30!     and consequently
31!            h=int_from_0_to_T cp*dT
32!     For the secant value one has:
33!            h=cp_secant*T
34!
35      integer imat,ntmat_,id,nshcon(*),four,i
36!
37      real*8 t1l,shcon(0:3,ntmat_,*),cp,physcon(*)
38!
39      four=4
40!
41!     calculating the tangent specific heat
42!
43      call ident2(shcon(0,1,imat),t1l,nshcon(imat),four,id)
44      if(nshcon(imat).eq.0) then
45         continue
46      elseif(nshcon(imat).eq.1) then
47         cp=shcon(1,1,imat)
48      elseif(id.eq.0) then
49         cp=shcon(1,1,imat)
50      elseif(id.eq.nshcon(imat)) then
51         cp=(shcon(0,1,imat)-physcon(1))*shcon(1,1,imat)
52         do i=2,nshcon(imat)
53            cp=cp+(shcon(0,i,imat)-shcon(0,i-1,imat))*
54     &              (shcon(1,i,imat)+shcon(1,i-1,imat))/2.d0
55         enddo
56         cp=cp+(t1l-shcon(0,nshcon(imat),imat))*
57     &           (shcon(1,nshcon(imat),imat))
58         cp=cp/(t1l-physcon(1))
59      else
60         cp=shcon(1,id,imat)+
61     &        (shcon(1,id+1,imat)-shcon(1,id,imat))*
62     &        (t1l-shcon(0,id,imat))/
63     &        (shcon(0,id+1,imat)-shcon(0,id,imat))
64         cp=(t1l-shcon(0,id,imat))*(cp+shcon(1,id,imat))/2.d0
65         do i=2,id
66            cp=cp+(shcon(0,i,imat)-shcon(0,i-1,imat))*
67     &              (shcon(1,i,imat)+shcon(1,i-1,imat))/2.d0
68         enddo
69         cp=cp+(shcon(0,1,imat)-physcon(1))*shcon(1,1,imat)
70         cp=cp/(t1l-physcon(1))
71      endif
72!
73      return
74      end
75
76
77
78
79
80
81
82