1!-------------------------------------------------------------------------------
2
3! This file is part of Code_Saturne, a general-purpose CFD tool.
4!
5! Copyright (C) 1998-2021 EDF S.A.
6!
7! This program is free software; you can redistribute it and/or modify it under
8! the terms of the GNU General Public License as published by the Free Software
9! Foundation; either version 2 of the License, or (at your option) any later
10! version.
11!
12! This program is distributed in the hope that it will be useful, but WITHOUT
13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
15! details.
16!
17! You should have received a copy of the GNU General Public License along with
18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19! Street, Fifth Floor, Boston, MA 02110-1301, USA.
20
21!-------------------------------------------------------------------------------
22!> \file rayive.f90
23!> \brief 1D Radiative scheme - IR H20 and dimere absorption
24
25!> \brief Compute in infrared spectral domain the emeissivity of H2O and dimere
26!-------------------------------------------------------------------------------
27! Arguments
28!______________________________________________________________________________.
29!  mode           name          role
30!______________________________________________________________________________!
31!> \param[out]      tauv      transmission function for water vapor and dimers
32!> \param[out]      dtauv     d(tauv)/dz
33!> \param[in]       qqqq      optical dpeth for water vapor (z,z')
34!> \param[in]       xqx       effective concentration absorption by water vapor
35!> \param[in]       qqqqc     idem qqqq for dimers
36!> \param[in]       xqc       idem xqx for dimers
37!> \param[in]       ro        air density
38!-------------------------------------------------------------------------------
39subroutine rayive (tauv,dtauv,qqqq,xqx,qqqqc,xqc,ro)
40
41implicit none
42
43!===============================================================================
44!... declaration des variables externes
45
46double precision tauv,dtauv,qqqq,xqx,qqqqc,xqc,ro
47
48!... declaration des variables internes
49
50double precision u,us,xu,xus,za,dza
51double precision a0,a1,a2,a3,a4,b0,b1,b2,b3,b4
52double precision as0,as1,as2,bs0,bs1,bs2
53double precision abs0,dabs0
54double precision n1,dn1,d1,dd1,t1,dt1,n2,dn2,d2,dd2,t2,dt2
55
56data a0,a1,a2,a3,a4/7.76192d-7,1.33836d-3,1.66649d-1,2.17686      &
57,2.69020/
58data b0,b1,b2,b3,b4/7.79097d-7,1.36832d-3,1.79601d-1,2.70573      &
59,5.15119/
60data as0,as1,as2/1.5075d-2,-3.6185d-2,1.9245d-2/
61data bs0,bs1,bs2/1.5075d-2,1.9547d-1,7.5271d-1/
62
63!===============================================================================
64
65u = qqqq/10.d0
66us = qqqqc/10.d0
67xu = xqx/10.d0
68xus = xqc/10.d0
69
70!  emissivite de la vapeur d'eau
71
72if(u.ge.0.01d0) then
73  za = 0.24d0*log10(u+0.01d0) + 0.622d0
74  dza = 0.24d0/(u + 0.01d0)/log(10.d0)
75else
76  za = 0.846d0*(u + 3.59d-5)**0.243d0 - 6.9d-2
77  dza = 0.846d0*0.243d0*(u + 3.59d-5)**(0.243d0 - 1.d0)
78endif
79
80!  emissivite des dimeres
81
82n1 = a0 + u*(a1 + u*(a2 + u*(a3 + u*a4)))
83dn1 = a1 + u*(2.d0*a2 + u*(3.d0*a3 + u*4.d0*a4))
84d1 = b0 + u*(b1 + u*(b2 + u*(b3 + u*(b4 + u))))
85dd1 = b1 + u*(2.d0*b2 + u*(3.d0*b3 + u*(4.d0*b4 + u*5.d0)))
86t1 = n1/d1
87
88dt1 = dn1/d1 - n1*dd1/d1/d1
89
90if(us.le.0.5d0) then
91  n2 = as0 + us*(as1 + us*as2)
92  dn2 = as1 + us*2.d0*as2
93  d2 = bs0 + us*(bs1 + us*(bs2 + us))
94  dd2 = bs1 + us*(2.d0*bs2 + 3.d0*us)
95  t2 = n2/d2
96
97  dt2 = dn2/d2 - n2*dd2/d2/d2
98else
99  t2 = 0.d0
100  dt2 = 0.d0
101endif
102
103!  transmission totale
104
105abs0 = za + 0.4614d0*t1*(1.d0 - t2)
106tauv = 1.d0 - abs0
107dabs0 = (dza*xu + 0.4614d0*(dt1*xu*(1.d0 - t2) - t1*dt2*xus))*ro
108dtauv = -dabs0
109
110return
111end subroutine rayive
112