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!cd  incompressible for ASME nozzles eq 4a 4b
20!
21!     author: Yannick Muller
22!
23      subroutine pk_cdi_noz(reynolds,cdi_noz)
24!
25      implicit none
26!
27      real*8 reynolds,cdi_noz,ln_reynolds,cdi_noz_lr,
28     &     cdi_noz_hr,e,reynolds_cor
29!
30      if (reynolds.lt.40000d0) then
31!
32! formerly pk_cdi_noz_lr : for low Reynolds nsumber
33!
34         if (reynolds.eq.0d0) then
35            reynolds_cor=1.d0
36         else
37            reynolds_cor=reynolds
38         endif
39         e=2.718281828459045d0
40         ln_reynolds=log(reynolds_cor)/log(e)
41!
42         cdi_noz_lr=0.19436d0+0.152884d0*ln_reynolds
43     &        -0.0097785d0*ln_reynolds**2d0+0.00020903d0
44     &        *ln_reynolds**3d0
45!
46         cdi_noz=cdi_noz_lr
47!
48      elseif (reynolds.lt.50000d0) then
49!
50         if (reynolds.eq.0) then
51            reynolds_cor=1
52         else
53            reynolds_cor=reynolds
54         endif
55!
56         e=2.718281828459045d0
57         ln_reynolds=log(reynolds_cor)/log(e)
58!
59         cdi_noz_lr=0.19436d0+0.152884d0*ln_reynolds
60     &        -0.0097785d0*ln_reynolds**2+0.00020903d0
61     &        *ln_reynolds**3d0
62!
63         cdi_noz_hr=0.9975d0-0.00653d0*dsqrt(1000000d0/50000d0)
64
65!     linear interpolation in order to achieve continuity
66!
67         cdi_noz=cdi_noz_lr+(cdi_noz_hr-cdi_noz_lr)
68     &        *(reynolds-40000d0)/(50000d0-40000d0)
69      else
70!
71!     formerly pk_cdi_noz_hr for high Reynolds numbers
72!
73         cdi_noz=0.9975d0-0.00653d0*dsqrt(1000000d0/reynolds)
74      endif
75
76      return
77      end
78