Line data Source code
1 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 : ! Copyright (c) 2015, Regents of the University of Colorado
3 : ! All rights reserved.
4 : !
5 : ! Redistribution and use in source and binary forms, with or without modification, are
6 : ! permitted provided that the following conditions are met:
7 : !
8 : ! 1. Redistributions of source code must retain the above copyright notice, this list of
9 : ! conditions and the following disclaimer.
10 : !
11 : ! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 : ! of conditions and the following disclaimer in the documentation and/or other
13 : ! materials provided with the distribution.
14 : !
15 : ! 3. Neither the name of the copyright holder nor the names of its contributors may be
16 : ! used to endorse or promote products derived from this software without specific prior
17 : ! written permission.
18 : !
19 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 : ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 : ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24 : ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 : ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 : ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 : ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 : !
29 : ! History:
30 : ! July 2006: John Haynes - Initial version
31 : ! May 2015: Dustin Swales - Modified for COSPv2.0
32 : !
33 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 : module optics_lib
35 : USE COSP_KINDS, ONLY: wp
36 : use mod_cosp_error, ONLY: errorMessage
37 : implicit none
38 :
39 : contains
40 :
41 : ! ##############################################################################
42 : ! subroutine M_WAT
43 : ! ##############################################################################
44 14316743 : subroutine m_wat(freq, tk, n_r, n_i)
45 : ! ############################################################################
46 : !
47 : ! Purpose:
48 : ! compute complex index of refraction of liquid water
49 : !
50 : ! Inputs:
51 : ! [freq] frequency (GHz)
52 : ! [tk] temperature (K)
53 : !
54 : ! Outputs:
55 : ! [n_r] real part index of refraction
56 : ! [n_i] imaginary part index of refraction
57 : !
58 : ! Reference:
59 : ! Based on the work of Ray (1972)
60 : !
61 : ! Coded:
62 : ! 03/22/05 John Haynes (haynes@atmos.colostate.edu)
63 : ! ############################################################################
64 :
65 : ! INPUTS
66 : real(wp), intent(in) :: &
67 : freq, & ! Frequency (GHz)
68 : tk ! Temperature (K)
69 :
70 : ! OUTPUTS
71 : real(wp), intent(out) :: &
72 : n_r, & ! Real part of index of refraction
73 : n_i ! Imaginary part of index of refraction
74 :
75 : ! Internal variables
76 : real(wp) :: ld,es,ei,a,ls,sg,tm1,cos1,sin1,e_r,e_i,pi,tc
77 : complex(wp) :: e_comp, sq
78 :
79 14316743 : tc = tk - 273.15_wp
80 :
81 14316743 : ld = 100._wp*2.99792458E8_wp/(freq*1E9_wp)
82 : es = 78.54_wp*(1-(4.579E-3_wp*(tc-25._wp)+1.19E-5_wp*(tc-25._wp)**2 &
83 14316743 : -2.8E-8_wp*(tc-25._wp)**3))
84 14316743 : ei = 5.27137_wp+0.021647_wp*tc-0.00131198_wp*tc**2
85 14316743 : a = -(16.8129_wp/(tc+273._wp))+0.0609265_wp
86 14316743 : ls = 0.00033836_wp*exp(2513.98_wp/(tc+273._wp))
87 14316743 : sg = 12.5664E8_wp
88 :
89 14316743 : tm1 = (ls/ld)**(1-a)
90 14316743 : pi = acos(-1._wp)
91 14316743 : cos1 = cos(0.5_wp*a*pi)
92 14316743 : sin1 = sin(0.5_wp*a*pi)
93 :
94 14316743 : e_r = ei + (((es-ei)*(1.+tm1*sin1))/(1._wp+2*tm1*sin1+tm1**2))
95 : e_i = (((es-ei)*tm1*cos1)/(1._wp+2*tm1*sin1+tm1**2)) &
96 14316743 : +((sg*ld)/1.885E11_wp)
97 :
98 : !ds e_comp = cmplx(e_r,e_i,Kind=Kind(0d0))
99 14316743 : e_comp = cmplx(e_r,e_i,Kind=wp)
100 14316743 : sq = sqrt(e_comp)
101 :
102 14316743 : n_r = real(sq)
103 14316743 : n_i = aimag(sq)
104 :
105 14316743 : return
106 : end subroutine m_wat
107 :
108 : ! ############################################################################
109 : ! subroutine M_ICE
110 : ! ############################################################################
111 13633725 : subroutine m_ice(freq,t,n_r,n_i)
112 : ! ##########################################################################
113 : !
114 : ! Purpose:
115 : ! compute complex index of refraction of ice
116 : !
117 : ! Inputs:
118 : ! [freq] frequency (GHz)
119 : ! [t] temperature (K)
120 : !
121 : ! Outputs:
122 : ! [n_r] real part index of refraction
123 : ! [n_i] imaginary part index of refraction
124 : !
125 : ! Reference:
126 : ! Fortran 90 port from IDL of REFICE by Stephen G. Warren
127 : !
128 : ! Modified:
129 : ! 05/31/05 John Haynes (haynes@atmos.colostate.edu)
130 : ! ##########################################################################
131 :
132 : ! INPUTS
133 : real(wp), intent(in) :: &
134 : freq, & ! Frequency (GHz)
135 : t ! Temperature (K)
136 :
137 : ! OUTPUTS
138 : real(wp), intent(out) :: &
139 : n_r, & ! Real part of index of refraction
140 : n_i ! Imaginary part of index of refraction
141 :
142 : ! Internal variables
143 : integer :: i,lt1,lt2
144 : real(wp) :: alam,pi,t1,t2, &
145 : x,x1,x2,y,y1,y2,ylo,yhi,tk
146 :
147 :
148 : ! Parameters:
149 : integer,parameter :: &
150 : nwl = 468, & !
151 : nwlt = 62 !
152 : real(wp),parameter,dimension(4) :: &
153 : temref = [272.16,268.16,253.16,213.16]
154 : real(wp),parameter :: & !
155 : wlmin = 0.045, & !
156 : wlmax = 8.6e6, & !
157 : cutice = 167.0
158 : real(wp),parameter,dimension(nwlt) :: &
159 : wlt = &
160 : [0.1670e+03, 0.1778e+03, 0.1884e+03, 0.1995e+03, 0.2113e+03, 0.2239e+03, &
161 : 0.2371e+03, 0.2512e+03, 0.2661e+03, 0.2818e+03, 0.2985e+03, 0.3162e+03, &
162 : 0.3548e+03, 0.3981e+03, 0.4467e+03, 0.5012e+03, 0.5623e+03, 0.6310e+03, &
163 : 0.7943e+03, 0.1000e+04, 0.1259e+04, 0.2500e+04, 0.5000e+04, 0.1000e+05, &
164 : 0.2000e+05, 0.3200e+05, 0.3500e+05, 0.4000e+05, 0.4500e+05, 0.5000e+05, &
165 : 0.6000e+05, 0.7000e+05, 0.9000e+05, 0.1110e+06, 0.1200e+06, 0.1300e+06, &
166 : 0.1400e+06, 0.1500e+06, 0.1600e+06, 0.1700e+06, 0.1800e+06, 0.2000e+06, &
167 : 0.2500e+06, 0.2900e+06, 0.3200e+06, 0.3500e+06, 0.3800e+06, 0.4000e+06, &
168 : 0.4500e+06, 0.5000e+06, 0.6000e+06, 0.6400e+06, 0.6800e+06, 0.7200e+06, &
169 : 0.7600e+06, 0.8000e+06, 0.8400e+06, 0.9000e+06, 0.1000e+07, 0.2000e+07, &
170 : 0.5000e+07,0.8600e+07]
171 : real(wp),parameter,dimension(nwl) :: &
172 : tabim = &
173 : [0.1640e+00, 0.1730e+00, 0.1830e+00, 0.1950e+00, 0.2080e+00, 0.2230e+00, &
174 : 0.2400e+00, 0.2500e+00, 0.2590e+00, 0.2680e+00, 0.2790e+00, 0.2970e+00, &
175 : 0.3190e+00, 0.3400e+00, 0.3660e+00, 0.3920e+00, 0.4160e+00, 0.4400e+00, &
176 : 0.4640e+00, 0.4920e+00, 0.5170e+00, 0.5280e+00, 0.5330e+00, 0.5340e+00, &
177 : 0.5310e+00, 0.5240e+00, 0.5100e+00, 0.5000e+00, 0.4990e+00, 0.4680e+00, &
178 : 0.3800e+00, 0.3600e+00, 0.3390e+00, 0.3180e+00, 0.2910e+00, 0.2510e+00, &
179 : 0.2440e+00, 0.2390e+00, 0.2390e+00, 0.2440e+00, 0.2470e+00, 0.2240e+00, &
180 : 0.1950e+00, 0.1740e+00, 0.1720e+00, 0.1800e+00, 0.1940e+00, 0.2130e+00, &
181 : 0.2430e+00, 0.2710e+00, 0.2890e+00, 0.3340e+00, 0.3440e+00, 0.3820e+00, &
182 : 0.4010e+00, 0.4065e+00, 0.4050e+00, 0.3890e+00, 0.3770e+00, 0.3450e+00, &
183 : 0.3320e+00, 0.3150e+00, 0.2980e+00, 0.2740e+00, 0.2280e+00, 0.1980e+00, &
184 : 0.1720e+00, 0.1560e+00, 0.1100e+00, 0.8300e-01, 0.5800e-01, 0.2200e-01, &
185 : 0.1000e-01, 0.3000e-02, 0.1000e-02, 0.3000e-03, 0.1000e-03, 0.3000e-04, &
186 : 0.1000e-04, 0.3000e-05, 0.1000e-05, 0.7000e-06, 0.4000e-06, 0.2000e-06, &
187 : 0.1000e-06, 0.6377e-07, 0.3750e-07, 0.2800e-07, 0.2400e-07, 0.2200e-07, &
188 : 0.1900e-07, 0.1750e-07, 0.1640e-07, 0.1590e-07, 0.1325e-07, 0.8623e-08, &
189 : 0.5504e-08, 0.3765e-08, 0.2710e-08, 0.2510e-08, 0.2260e-08, 0.2080e-08, &
190 : 0.1910e-08, 0.1540e-08, 0.1530e-08, 0.1550e-08, 0.1640e-08, 0.1780e-08, &
191 : 0.1910e-08, 0.2140e-08, 0.2260e-08, 0.2540e-08, 0.2930e-08, 0.3110e-08, &
192 : 0.3290e-08, 0.3520e-08, 0.4040e-08, 0.4880e-08, 0.5730e-08, 0.6890e-08, &
193 : 0.8580e-08, 0.1040e-07, 0.1220e-07, 0.1430e-07, 0.1660e-07, 0.1890e-07, &
194 : 0.2090e-07, 0.2400e-07, 0.2900e-07, 0.3440e-07, 0.4030e-07, 0.4300e-07, &
195 : 0.4920e-07, 0.5870e-07, 0.7080e-07, 0.8580e-07, 0.1020e-06, 0.1180e-06, &
196 : 0.1340e-06, 0.1400e-06, 0.1430e-06, 0.1450e-06, 0.1510e-06, 0.1830e-06, &
197 : 0.2150e-06, 0.2650e-06, 0.3350e-06, 0.3920e-06, 0.4200e-06, 0.4440e-06, &
198 : 0.4740e-06, 0.5110e-06, 0.5530e-06, 0.6020e-06, 0.7550e-06, 0.9260e-06, &
199 : 0.1120e-05, 0.1330e-05, 0.1620e-05, 0.2000e-05, 0.2250e-05, 0.2330e-05, &
200 : 0.2330e-05, 0.2170e-05, 0.1960e-05, 0.1810e-05, 0.1740e-05, 0.1730e-05, &
201 : 0.1700e-05, 0.1760e-05, 0.1820e-05, 0.2040e-05, 0.2250e-05, 0.2290e-05, &
202 : 0.3040e-05, 0.3840e-05, 0.4770e-05, 0.5760e-05, 0.6710e-05, 0.8660e-05, &
203 : 0.1020e-04, 0.1130e-04, 0.1220e-04, 0.1290e-04, 0.1320e-04, 0.1350e-04, &
204 : 0.1330e-04, 0.1320e-04, 0.1320e-04, 0.1310e-04, 0.1320e-04, 0.1320e-04, &
205 : 0.1340e-04, 0.1390e-04, 0.1420e-04, 0.1480e-04, 0.1580e-04, 0.1740e-04, &
206 : 0.1980e-04, 0.2500e-04, 0.5400e-04, 0.1040e-03, 0.2030e-03, 0.2708e-03, &
207 : 0.3511e-03, 0.4299e-03, 0.5181e-03, 0.5855e-03, 0.5899e-03, 0.5635e-03, &
208 : 0.5480e-03, 0.5266e-03, 0.4394e-03, 0.3701e-03, 0.3372e-03, 0.2410e-03, &
209 : 0.1890e-03, 0.1660e-03, 0.1450e-03, 0.1280e-03, 0.1030e-03, 0.8600e-04, &
210 : 0.8220e-04, 0.8030e-04, 0.8500e-04, 0.9900e-04, 0.1500e-03, 0.2950e-03, &
211 : 0.4687e-03, 0.7615e-03, 0.1010e-02, 0.1313e-02, 0.1539e-02, 0.1588e-02, &
212 : 0.1540e-02, 0.1412e-02, 0.1244e-02, 0.1068e-02, 0.8414e-03, 0.5650e-03, &
213 : 0.4320e-03, 0.3500e-03, 0.2870e-03, 0.2210e-03, 0.2030e-03, 0.2010e-03, &
214 : 0.2030e-03, 0.2140e-03, 0.2320e-03, 0.2890e-03, 0.3810e-03, 0.4620e-03, &
215 : 0.5480e-03, 0.6180e-03, 0.6800e-03, 0.7300e-03, 0.7820e-03, 0.8480e-03, &
216 : 0.9250e-03, 0.9200e-03, 0.8920e-03, 0.8700e-03, 0.8900e-03, 0.9300e-03, &
217 : 0.1010e-02, 0.1350e-02, 0.3420e-02, 0.7920e-02, 0.2000e-01, 0.3800e-01, &
218 : 0.5200e-01, 0.6800e-01, 0.9230e-01, 0.1270e+00, 0.1690e+00, 0.2210e+00, &
219 : 0.2760e+00, 0.3120e+00, 0.3470e+00, 0.3880e+00, 0.4380e+00, 0.4930e+00, &
220 : 0.5540e+00, 0.6120e+00, 0.6250e+00, 0.5930e+00, 0.5390e+00, 0.4910e+00, &
221 : 0.4380e+00, 0.3720e+00, 0.3000e+00, 0.2380e+00, 0.1930e+00, 0.1580e+00, &
222 : 0.1210e+00, 0.1030e+00, 0.8360e-01, 0.6680e-01, 0.5400e-01, 0.4220e-01, &
223 : 0.3420e-01, 0.2740e-01, 0.2200e-01, 0.1860e-01, 0.1520e-01, 0.1260e-01, &
224 : 0.1060e-01, 0.8020e-02, 0.6850e-02, 0.6600e-02, 0.6960e-02, 0.9160e-02, &
225 : 0.1110e-01, 0.1450e-01, 0.2000e-01, 0.2300e-01, 0.2600e-01, 0.2900e-01, &
226 : 0.2930e-01, 0.3000e-01, 0.2850e-01, 0.1730e-01, 0.1290e-01, 0.1200e-01, &
227 : 0.1250e-01, 0.1340e-01, 0.1400e-01, 0.1750e-01, 0.2400e-01, 0.3500e-01, &
228 : 0.3800e-01, 0.4200e-01, 0.4600e-01, 0.5200e-01, 0.5700e-01, 0.6900e-01, &
229 : 0.7000e-01, 0.6700e-01, 0.6500e-01, 0.6400e-01, 0.6200e-01, 0.5900e-01, &
230 : 0.5700e-01, 0.5600e-01, 0.5500e-01, 0.5700e-01, 0.5800e-01, 0.5700e-01, &
231 : 0.5500e-01, 0.5500e-01, 0.5400e-01, 0.5200e-01, 0.5200e-01, 0.5200e-01, &
232 : 0.5200e-01, 0.5000e-01, 0.4700e-01, 0.4300e-01, 0.3900e-01, 0.3700e-01, &
233 : 0.3900e-01, 0.4000e-01, 0.4200e-01, 0.4400e-01, 0.4500e-01, 0.4600e-01, &
234 : 0.4700e-01, 0.5100e-01, 0.6500e-01, 0.7500e-01, 0.8800e-01, 0.1080e+00, &
235 : 0.1340e+00, 0.1680e+00, 0.2040e+00, 0.2480e+00, 0.2800e+00, 0.3410e+00, &
236 : 0.3790e+00, 0.4090e+00, 0.4220e+00, 0.4220e+00, 0.4030e+00, 0.3890e+00, &
237 : 0.3740e+00, 0.3540e+00, 0.3350e+00, 0.3150e+00, 0.2940e+00, 0.2710e+00, &
238 : 0.2460e+00, 0.1980e+00, 0.1640e+00, 0.1520e+00, 0.1420e+00, 0.1280e+00, &
239 : 0.1250e+00, 0.1230e+00, 0.1160e+00, 0.1070e+00, 0.7900e-01, 0.7200e-01, &
240 : 0.7600e-01, 0.7500e-01, 0.6700e-01, 0.5500e-01, 0.4500e-01, 0.2900e-01, &
241 : 0.2750e-01, 0.2700e-01, 0.2730e-01, 0.2890e-01, 0.3000e-01, 0.3400e-01, &
242 : 0.5300e-01, 0.7550e-01, 0.1060e+00, 0.1350e+00, 0.1761e+00, 0.2229e+00, &
243 : 0.2746e+00, 0.3280e+00, 0.3906e+00, 0.4642e+00, 0.5247e+00, 0.5731e+00, &
244 : 0.6362e+00, 0.6839e+00, 0.7091e+00, 0.6790e+00, 0.6250e+00, 0.5654e+00, &
245 : 0.5433e+00, 0.5292e+00, 0.5070e+00, 0.4883e+00, 0.4707e+00, 0.4203e+00, &
246 : 0.3771e+00, 0.3376e+00, 0.3056e+00, 0.2835e+00, 0.3170e+00, 0.3517e+00, &
247 : 0.3902e+00, 0.4509e+00, 0.4671e+00, 0.4779e+00, 0.4890e+00, 0.4899e+00, &
248 : 0.4873e+00, 0.4766e+00, 0.4508e+00, 0.4193e+00, 0.3880e+00, 0.3433e+00, &
249 : 0.3118e+00, 0.2935e+00, 0.2350e+00, 0.1981e+00, 0.1865e+00, 0.1771e+00, &
250 : 0.1620e+00, 0.1490e+00, 0.1390e+00, 0.1200e+00, 0.9620e-01, 0.8300e-01]
251 : real(wp),parameter,dimension(nwl) :: &
252 : wl = &
253 : [0.4430e-01, 0.4510e-01, 0.4590e-01, 0.4680e-01, 0.4770e-01, 0.4860e-01, &
254 : 0.4960e-01, 0.5060e-01, 0.5170e-01, 0.5280e-01, 0.5390e-01, 0.5510e-01, &
255 : 0.5640e-01, 0.5770e-01, 0.5900e-01, 0.6050e-01, 0.6200e-01, 0.6360e-01, &
256 : 0.6530e-01, 0.6700e-01, 0.6890e-01, 0.7080e-01, 0.7290e-01, 0.7380e-01, &
257 : 0.7510e-01, 0.7750e-01, 0.8000e-01, 0.8270e-01, 0.8550e-01, 0.8860e-01, &
258 : 0.9180e-01, 0.9300e-01, 0.9540e-01, 0.9920e-01, 0.1033e+00, 0.1078e+00, &
259 : 0.1100e+00, 0.1127e+00, 0.1140e+00, 0.1181e+00, 0.1210e+00, 0.1240e+00, &
260 : 0.1272e+00, 0.1295e+00, 0.1305e+00, 0.1319e+00, 0.1333e+00, 0.1348e+00, &
261 : 0.1362e+00, 0.1370e+00, 0.1378e+00, 0.1387e+00, 0.1393e+00, 0.1409e+00, &
262 : 0.1425e+00, 0.1435e+00, 0.1442e+00, 0.1450e+00, 0.1459e+00, 0.1468e+00, &
263 : 0.1476e+00, 0.1480e+00, 0.1485e+00, 0.1494e+00, 0.1512e+00, 0.1531e+00, &
264 : 0.1540e+00, 0.1550e+00, 0.1569e+00, 0.1580e+00, 0.1589e+00, 0.1610e+00, &
265 : 0.1625e+00, 0.1648e+00, 0.1669e+00, 0.1692e+00, 0.1713e+00, 0.1737e+00, &
266 : 0.1757e+00, 0.1779e+00, 0.1802e+00, 0.1809e+00, 0.1821e+00, 0.1833e+00, &
267 : 0.1843e+00, 0.1850e+00, 0.1860e+00, 0.1870e+00, 0.1880e+00, 0.1890e+00, &
268 : 0.1900e+00, 0.1910e+00, 0.1930e+00, 0.1950e+00, 0.2100e+00, 0.2500e+00, &
269 : 0.3000e+00, 0.3500e+00, 0.4000e+00, 0.4100e+00, 0.4200e+00, 0.4300e+00, &
270 : 0.4400e+00, 0.4500e+00, 0.4600e+00, 0.4700e+00, 0.4800e+00, 0.4900e+00, &
271 : 0.5000e+00, 0.5100e+00, 0.5200e+00, 0.5300e+00, 0.5400e+00, 0.5500e+00, &
272 : 0.5600e+00, 0.5700e+00, 0.5800e+00, 0.5900e+00, 0.6000e+00, 0.6100e+00, &
273 : 0.6200e+00, 0.6300e+00, 0.6400e+00, 0.6500e+00, 0.6600e+00, 0.6700e+00, &
274 : 0.6800e+00, 0.6900e+00, 0.7000e+00, 0.7100e+00, 0.7200e+00, 0.7300e+00, &
275 : 0.7400e+00, 0.7500e+00, 0.7600e+00, 0.7700e+00, 0.7800e+00, 0.7900e+00, &
276 : 0.8000e+00, 0.8100e+00, 0.8200e+00, 0.8300e+00, 0.8400e+00, 0.8500e+00, &
277 : 0.8600e+00, 0.8700e+00, 0.8800e+00, 0.8900e+00, 0.9000e+00, 0.9100e+00, &
278 : 0.9200e+00, 0.9300e+00, 0.9400e+00, 0.9500e+00, 0.9600e+00, 0.9700e+00, &
279 : 0.9800e+00, 0.9900e+00, 0.1000e+01, 0.1010e+01, 0.1020e+01, 0.1030e+01, &
280 : 0.1040e+01, 0.1050e+01, 0.1060e+01, 0.1070e+01, 0.1080e+01, 0.1090e+01, &
281 : 0.1100e+01, 0.1110e+01, 0.1120e+01, 0.1130e+01, 0.1140e+01, 0.1150e+01, &
282 : 0.1160e+01, 0.1170e+01, 0.1180e+01, 0.1190e+01, 0.1200e+01, 0.1210e+01, &
283 : 0.1220e+01, 0.1230e+01, 0.1240e+01, 0.1250e+01, 0.1260e+01, 0.1270e+01, &
284 : 0.1280e+01, 0.1290e+01, 0.1300e+01, 0.1310e+01, 0.1320e+01, 0.1330e+01, &
285 : 0.1340e+01, 0.1350e+01, 0.1360e+01, 0.1370e+01, 0.1380e+01, 0.1390e+01, &
286 : 0.1400e+01, 0.1410e+01, 0.1420e+01, 0.1430e+01, 0.1440e+01, 0.1449e+01, &
287 : 0.1460e+01, 0.1471e+01, 0.1481e+01, 0.1493e+01, 0.1504e+01, 0.1515e+01, &
288 : 0.1527e+01, 0.1538e+01, 0.1563e+01, 0.1587e+01, 0.1613e+01, 0.1650e+01, &
289 : 0.1680e+01, 0.1700e+01, 0.1730e+01, 0.1760e+01, 0.1800e+01, 0.1830e+01, &
290 : 0.1840e+01, 0.1850e+01, 0.1855e+01, 0.1860e+01, 0.1870e+01, 0.1890e+01, &
291 : 0.1905e+01, 0.1923e+01, 0.1942e+01, 0.1961e+01, 0.1980e+01, 0.2000e+01, &
292 : 0.2020e+01, 0.2041e+01, 0.2062e+01, 0.2083e+01, 0.2105e+01, 0.2130e+01, &
293 : 0.2150e+01, 0.2170e+01, 0.2190e+01, 0.2220e+01, 0.2240e+01, 0.2245e+01, &
294 : 0.2250e+01, 0.2260e+01, 0.2270e+01, 0.2290e+01, 0.2310e+01, 0.2330e+01, &
295 : 0.2350e+01, 0.2370e+01, 0.2390e+01, 0.2410e+01, 0.2430e+01, 0.2460e+01, &
296 : 0.2500e+01, 0.2520e+01, 0.2550e+01, 0.2565e+01, 0.2580e+01, 0.2590e+01, &
297 : 0.2600e+01, 0.2620e+01, 0.2675e+01, 0.2725e+01, 0.2778e+01, 0.2817e+01, &
298 : 0.2833e+01, 0.2849e+01, 0.2865e+01, 0.2882e+01, 0.2899e+01, 0.2915e+01, &
299 : 0.2933e+01, 0.2950e+01, 0.2967e+01, 0.2985e+01, 0.3003e+01, 0.3021e+01, &
300 : 0.3040e+01, 0.3058e+01, 0.3077e+01, 0.3096e+01, 0.3115e+01, 0.3135e+01, &
301 : 0.3155e+01, 0.3175e+01, 0.3195e+01, 0.3215e+01, 0.3236e+01, 0.3257e+01, &
302 : 0.3279e+01, 0.3300e+01, 0.3322e+01, 0.3345e+01, 0.3367e+01, 0.3390e+01, &
303 : 0.3413e+01, 0.3436e+01, 0.3460e+01, 0.3484e+01, 0.3509e+01, 0.3534e+01, &
304 : 0.3559e+01, 0.3624e+01, 0.3732e+01, 0.3775e+01, 0.3847e+01, 0.3969e+01, &
305 : 0.4099e+01, 0.4239e+01, 0.4348e+01, 0.4387e+01, 0.4444e+01, 0.4505e+01, &
306 : 0.4547e+01, 0.4560e+01, 0.4580e+01, 0.4719e+01, 0.4904e+01, 0.5000e+01, &
307 : 0.5100e+01, 0.5200e+01, 0.5263e+01, 0.5400e+01, 0.5556e+01, 0.5714e+01, &
308 : 0.5747e+01, 0.5780e+01, 0.5814e+01, 0.5848e+01, 0.5882e+01, 0.6061e+01, &
309 : 0.6135e+01, 0.6250e+01, 0.6289e+01, 0.6329e+01, 0.6369e+01, 0.6410e+01, &
310 : 0.6452e+01, 0.6494e+01, 0.6579e+01, 0.6667e+01, 0.6757e+01, 0.6897e+01, &
311 : 0.7042e+01, 0.7143e+01, 0.7246e+01, 0.7353e+01, 0.7463e+01, 0.7576e+01, &
312 : 0.7692e+01, 0.7812e+01, 0.7937e+01, 0.8065e+01, 0.8197e+01, 0.8333e+01, &
313 : 0.8475e+01, 0.8696e+01, 0.8929e+01, 0.9091e+01, 0.9259e+01, 0.9524e+01, &
314 : 0.9804e+01, 0.1000e+02, 0.1020e+02, 0.1031e+02, 0.1042e+02, 0.1053e+02, &
315 : 0.1064e+02, 0.1075e+02, 0.1087e+02, 0.1100e+02, 0.1111e+02, 0.1136e+02, &
316 : 0.1163e+02, 0.1190e+02, 0.1220e+02, 0.1250e+02, 0.1282e+02, 0.1299e+02, &
317 : 0.1316e+02, 0.1333e+02, 0.1351e+02, 0.1370e+02, 0.1389e+02, 0.1408e+02, &
318 : 0.1429e+02, 0.1471e+02, 0.1515e+02, 0.1538e+02, 0.1563e+02, 0.1613e+02, &
319 : 0.1639e+02, 0.1667e+02, 0.1695e+02, 0.1724e+02, 0.1818e+02, 0.1887e+02, &
320 : 0.1923e+02, 0.1961e+02, 0.2000e+02, 0.2041e+02, 0.2083e+02, 0.2222e+02, &
321 : 0.2260e+02, 0.2305e+02, 0.2360e+02, 0.2460e+02, 0.2500e+02, 0.2600e+02, &
322 : 0.2857e+02, 0.3100e+02, 0.3333e+02, 0.3448e+02, 0.3564e+02, 0.3700e+02, &
323 : 0.3824e+02, 0.3960e+02, 0.4114e+02, 0.4276e+02, 0.4358e+02, 0.4458e+02, &
324 : 0.4550e+02, 0.4615e+02, 0.4671e+02, 0.4736e+02, 0.4800e+02, 0.4878e+02, &
325 : 0.5003e+02, 0.5128e+02, 0.5275e+02, 0.5350e+02, 0.5424e+02, 0.5500e+02, &
326 : 0.5574e+02, 0.5640e+02, 0.5700e+02, 0.5746e+02, 0.5840e+02, 0.5929e+02, &
327 : 0.6000e+02, 0.6100e+02, 0.6125e+02, 0.6250e+02, 0.6378e+02, 0.6467e+02, &
328 : 0.6558e+02, 0.6655e+02, 0.6760e+02, 0.6900e+02, 0.7053e+02, 0.7300e+02, &
329 : 0.7500e+02, 0.7629e+02, 0.8000e+02, 0.8297e+02, 0.8500e+02, 0.8680e+02, &
330 : 0.9080e+02, 0.9517e+02, 0.1000e+03, 0.1200e+03, 0.1500e+03, 0.1670e+03]
331 : real(wp),parameter,dimension(nwlt,4) :: &
332 : tabimt = reshape(source= &
333 : (/0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4560e-01, 0.3790e-01, 0.3140e-01, &
334 : 0.2620e-01, 0.2240e-01, 0.1960e-01, 0.1760e-01, 0.1665e-01, 0.1620e-01, &
335 : 0.1550e-01, 0.1470e-01, 0.1390e-01, 0.1320e-01, 0.1250e-01, 0.1180e-01, &
336 : 0.1060e-01, 0.9540e-02, 0.8560e-02, 0.6210e-02, 0.4490e-02, 0.3240e-02, &
337 : 0.2340e-02, 0.1880e-02, 0.1740e-02, 0.1500e-02, 0.1320e-02, 0.1160e-02, &
338 : 0.8800e-03, 0.6950e-03, 0.4640e-03, 0.3400e-03, 0.3110e-03, 0.2940e-03, &
339 : 0.2790e-03, 0.2700e-03, 0.2640e-03, 0.2580e-03, 0.2520e-03, 0.2490e-03, &
340 : 0.2540e-03, 0.2640e-03, 0.2740e-03, 0.2890e-03, 0.3050e-03, 0.3150e-03, &
341 : 0.3460e-03, 0.3820e-03, 0.4620e-03, 0.5000e-03, 0.5500e-03, 0.5950e-03, &
342 : 0.6470e-03, 0.6920e-03, 0.7420e-03, 0.8200e-03, 0.9700e-03, 0.1950e-02, &
343 : 0.5780e-02, 0.9700e-02, 0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4560e-01, &
344 : 0.3790e-01, 0.3140e-01, 0.2620e-01, 0.2240e-01, 0.1960e-01, 0.1760e-01, &
345 : 0.1665e-01, 0.1600e-01, 0.1500e-01, 0.1400e-01, 0.1310e-01, 0.1230e-01, &
346 : 0.1150e-01, 0.1080e-01, 0.9460e-02, 0.8290e-02, 0.7270e-02, 0.4910e-02, &
347 : 0.3300e-02, 0.2220e-02, 0.1490e-02, 0.1140e-02, 0.1060e-02, 0.9480e-03, &
348 : 0.8500e-03, 0.7660e-03, 0.6300e-03, 0.5200e-03, 0.3840e-03, 0.2960e-03, &
349 : 0.2700e-03, 0.2520e-03, 0.2440e-03, 0.2360e-03, 0.2300e-03, 0.2280e-03, &
350 : 0.2250e-03, 0.2200e-03, 0.2160e-03, 0.2170e-03, 0.2200e-03, 0.2250e-03, &
351 : 0.2320e-03, 0.2390e-03, 0.2600e-03, 0.2860e-03, 0.3560e-03, 0.3830e-03, &
352 : 0.4150e-03, 0.4450e-03, 0.4760e-03, 0.5080e-03, 0.5400e-03, 0.5860e-03, &
353 : 0.6780e-03, 0.1280e-02, 0.3550e-02, 0.5600e-02, 0.8300e-01, 0.6900e-01, &
354 : 0.5700e-01, 0.4560e-01, 0.3790e-01, 0.3140e-01, 0.2620e-01, 0.2190e-01, &
355 : 0.1880e-01, 0.1660e-01, 0.1540e-01, 0.1470e-01, 0.1350e-01, 0.1250e-01, &
356 : 0.1150e-01, 0.1060e-01, 0.9770e-02, 0.9010e-02, 0.7660e-02, 0.6520e-02, &
357 : 0.5540e-02, 0.3420e-02, 0.2100e-02, 0.1290e-02, 0.7930e-03, 0.5700e-03, &
358 : 0.5350e-03, 0.4820e-03, 0.4380e-03, 0.4080e-03, 0.3500e-03, 0.3200e-03, &
359 : 0.2550e-03, 0.2120e-03, 0.2000e-03, 0.1860e-03, 0.1750e-03, 0.1660e-03, &
360 : 0.1560e-03, 0.1490e-03, 0.1440e-03, 0.1350e-03, 0.1210e-03, 0.1160e-03, &
361 : 0.1160e-03, 0.1170e-03, 0.1200e-03, 0.1230e-03, 0.1320e-03, 0.1440e-03, &
362 : 0.1680e-03, 0.1800e-03, 0.1900e-03, 0.2090e-03, 0.2160e-03, 0.2290e-03, &
363 : 0.2400e-03, 0.2600e-03, 0.2920e-03, 0.6100e-03, 0.1020e-02, 0.1810e-02, &
364 : 0.8300e-01, 0.6900e-01, 0.5700e-01, 0.4450e-01, 0.3550e-01, 0.2910e-01, &
365 : 0.2440e-01, 0.1970e-01, 0.1670e-01, 0.1400e-01, 0.1235e-01, 0.1080e-01, &
366 : 0.8900e-02, 0.7340e-02, 0.6400e-02, 0.5600e-02, 0.5000e-02, 0.4520e-02, &
367 : 0.3680e-02, 0.2990e-02, 0.2490e-02, 0.1550e-02, 0.9610e-03, 0.5950e-03, &
368 : 0.3690e-03, 0.2670e-03, 0.2510e-03, 0.2290e-03, 0.2110e-03, 0.1960e-03, &
369 : 0.1730e-03, 0.1550e-03, 0.1310e-03, 0.1130e-03, 0.1060e-03, 0.9900e-04, &
370 : 0.9300e-04, 0.8730e-04, 0.8300e-04, 0.7870e-04, 0.7500e-04, 0.6830e-04, &
371 : 0.5600e-04, 0.4960e-04, 0.4550e-04, 0.4210e-04, 0.3910e-04, 0.3760e-04, &
372 : 0.3400e-04, 0.3100e-04, 0.2640e-04, 0.2510e-04, 0.2430e-04, 0.2390e-04, &
373 : 0.2370e-04, 0.2380e-04, 0.2400e-04, 0.2460e-04, 0.2660e-04, 0.4450e-04, &
374 : 0.8700e-04, 0.1320e-03/),shape=(/nwlt,4/))
375 :
376 : real(wp),parameter,dimension(nwl) :: &
377 : tabre = &
378 : [0.83441, 0.83676, 0.83729, 0.83771, 0.83827, 0.84038, &
379 : 0.84719, 0.85522, 0.86047, 0.86248, 0.86157, 0.86093, &
380 : 0.86419, 0.86916, 0.87764, 0.89296, 0.91041, 0.93089, &
381 : 0.95373, 0.98188, 1.02334, 1.06735, 1.11197, 1.13134, &
382 : 1.15747, 1.20045, 1.23840, 1.27325, 1.32157, 1.38958, &
383 : 1.41644, 1.40906, 1.40063, 1.40169, 1.40934, 1.40221, &
384 : 1.39240, 1.38424, 1.38075, 1.38186, 1.39634, 1.40918, &
385 : 1.40256, 1.38013, 1.36303, 1.34144, 1.32377, 1.30605, &
386 : 1.29054, 1.28890, 1.28931, 1.30190, 1.32025, 1.36302, &
387 : 1.41872, 1.45834, 1.49028, 1.52128, 1.55376, 1.57782, &
388 : 1.59636, 1.60652, 1.61172, 1.61919, 1.62522, 1.63404, &
389 : 1.63689, 1.63833, 1.63720, 1.63233, 1.62222, 1.58269, &
390 : 1.55635, 1.52453, 1.50320, 1.48498, 1.47226, 1.45991, &
391 : 1.45115, 1.44272, 1.43498, 1.43280, 1.42924, 1.42602, &
392 : 1.42323, 1.42143, 1.41897, 1.41660, 1.41434, 1.41216, &
393 : 1.41006, 1.40805, 1.40423, 1.40067, 1.38004, 1.35085, &
394 : 1.33394, 1.32492, 1.31940, 1.31854, 1.31775, 1.31702, &
395 : 1.31633, 1.31569, 1.31509, 1.31452, 1.31399, 1.31349, &
396 : 1.31302, 1.31257, 1.31215, 1.31175, 1.31136, 1.31099, &
397 : 1.31064, 1.31031, 1.30999, 1.30968, 1.30938, 1.30909, &
398 : 1.30882, 1.30855, 1.30829, 1.30804, 1.30780, 1.30756, &
399 : 1.30733, 1.30710, 1.30688, 1.30667, 1.30646, 1.30625, &
400 : 1.30605, 1.30585, 1.30566, 1.30547, 1.30528, 1.30509, &
401 : 1.30491, 1.30473, 1.30455, 1.30437, 1.30419, 1.30402, &
402 : 1.30385, 1.30367, 1.30350, 1.30333, 1.30316, 1.30299, &
403 : 1.30283, 1.30266, 1.30249, 1.30232, 1.30216, 1.30199, &
404 : 1.30182, 1.30166, 1.30149, 1.30132, 1.30116, 1.30099, &
405 : 1.30082, 1.30065, 1.30048, 1.30031, 1.30014, 1.29997, &
406 : 1.29979, 1.29962, 1.29945, 1.29927, 1.29909, 1.29891, &
407 : 1.29873, 1.29855, 1.29837, 1.29818, 1.29800, 1.29781, &
408 : 1.29762, 1.29743, 1.29724, 1.29705, 1.29686, 1.29666, &
409 : 1.29646, 1.29626, 1.29605, 1.29584, 1.29563, 1.29542, &
410 : 1.29521, 1.29499, 1.29476, 1.29453, 1.29430, 1.29406, &
411 : 1.29381, 1.29355, 1.29327, 1.29299, 1.29272, 1.29252, &
412 : 1.29228, 1.29205, 1.29186, 1.29167, 1.29150, 1.29130, &
413 : 1.29106, 1.29083, 1.29025, 1.28962, 1.28891, 1.28784, &
414 : 1.28689, 1.28623, 1.28521, 1.28413, 1.28261, 1.28137, &
415 : 1.28093, 1.28047, 1.28022, 1.27998, 1.27948, 1.27849, &
416 : 1.27774, 1.27691, 1.27610, 1.27535, 1.27471, 1.27404, &
417 : 1.27329, 1.27240, 1.27139, 1.27029, 1.26901, 1.26736, &
418 : 1.26591, 1.26441, 1.26284, 1.26036, 1.25860, 1.25815, &
419 : 1.25768, 1.25675, 1.25579, 1.25383, 1.25179, 1.24967, &
420 : 1.24745, 1.24512, 1.24266, 1.24004, 1.23725, 1.23270, &
421 : 1.22583, 1.22198, 1.21548, 1.21184, 1.20790, 1.20507, &
422 : 1.20209, 1.19566, 1.17411, 1.14734, 1.10766, 1.06739, &
423 : 1.04762, 1.02650, 1.00357, 0.98197, 0.96503, 0.95962, &
424 : 0.97269, 0.99172, 1.00668, 1.02186, 1.04270, 1.07597, &
425 : 1.12954, 1.21267, 1.32509, 1.42599, 1.49656, 1.55095, &
426 : 1.59988, 1.63631, 1.65024, 1.64278, 1.62691, 1.61284, &
427 : 1.59245, 1.57329, 1.55770, 1.54129, 1.52654, 1.51139, &
428 : 1.49725, 1.48453, 1.47209, 1.46125, 1.45132, 1.44215, &
429 : 1.43366, 1.41553, 1.39417, 1.38732, 1.37735, 1.36448, &
430 : 1.35414, 1.34456, 1.33882, 1.33807, 1.33847, 1.34053, &
431 : 1.34287, 1.34418, 1.34634, 1.34422, 1.33453, 1.32897, &
432 : 1.32333, 1.31800, 1.31432, 1.30623, 1.29722, 1.28898, &
433 : 1.28730, 1.28603, 1.28509, 1.28535, 1.28813, 1.30156, &
434 : 1.30901, 1.31720, 1.31893, 1.32039, 1.32201, 1.32239, &
435 : 1.32149, 1.32036, 1.31814, 1.31705, 1.31807, 1.31953, &
436 : 1.31933, 1.31896, 1.31909, 1.31796, 1.31631, 1.31542, &
437 : 1.31540, 1.31552, 1.31455, 1.31193, 1.30677, 1.29934, &
438 : 1.29253, 1.28389, 1.27401, 1.26724, 1.25990, 1.24510, &
439 : 1.22241, 1.19913, 1.17150, 1.15528, 1.13700, 1.11808, &
440 : 1.10134, 1.09083, 1.08734, 1.09254, 1.10654, 1.14779, &
441 : 1.20202, 1.25825, 1.32305, 1.38574, 1.44478, 1.47170, &
442 : 1.49619, 1.51652, 1.53328, 1.54900, 1.56276, 1.57317, &
443 : 1.58028, 1.57918, 1.56672, 1.55869, 1.55081, 1.53807, &
444 : 1.53296, 1.53220, 1.53340, 1.53289, 1.51705, 1.50097, &
445 : 1.49681, 1.49928, 1.50153, 1.49856, 1.49053, 1.46070, &
446 : 1.45182, 1.44223, 1.43158, 1.41385, 1.40676, 1.38955, &
447 : 1.34894, 1.31039, 1.26420, 1.23656, 1.21663, 1.20233, &
448 : 1.19640, 1.19969, 1.20860, 1.22173, 1.24166, 1.28175, &
449 : 1.32784, 1.38657, 1.46486, 1.55323, 1.60379, 1.61877, &
450 : 1.62963, 1.65712, 1.69810, 1.72065, 1.74865, 1.76736, &
451 : 1.76476, 1.75011, 1.72327, 1.68490, 1.62398, 1.59596, &
452 : 1.58514, 1.59917, 1.61405, 1.66625, 1.70663, 1.73713, &
453 : 1.76860, 1.80343, 1.83296, 1.85682, 1.87411, 1.89110, &
454 : 1.89918, 1.90432, 1.90329, 1.88744, 1.87499, 1.86702, &
455 : 1.85361, 1.84250, 1.83225, 1.81914, 1.82268, 1.82961]
456 : real(wp),parameter,dimension(nwlt,4) :: &
457 : tabret = reshape( &
458 : source =(/1.82961, 1.83258, 1.83149, &
459 : 1.82748, 1.82224, 1.81718, 1.81204, 1.80704, 1.80250, &
460 : 1.79834, 1.79482, 1.79214, 1.78843, 1.78601, 1.78434, &
461 : 1.78322, 1.78248, 1.78201, 1.78170, 1.78160, 1.78190, &
462 : 1.78300, 1.78430, 1.78520, 1.78620, 1.78660, 1.78680, &
463 : 1.78690, 1.78700, 1.78700, 1.78710, 1.78710, 1.78720, &
464 : 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, &
465 : 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, &
466 : 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, &
467 : 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, 1.78720, &
468 : 1.78720, 1.78720, 1.78720, 1.78720, 1.78800, &
469 : 1.82961, 1.83258, 1.83149, 1.82748, &
470 : 1.82224, 1.81718, 1.81204, 1.80704, 1.80250, 1.79834, &
471 : 1.79482, 1.79214, 1.78843, 1.78601, 1.78434, 1.78322, &
472 : 1.78248, 1.78201, 1.78170, 1.78160, 1.78190, 1.78300, &
473 : 1.78430, 1.78520, 1.78610, 1.78630, 1.78640, 1.78650, &
474 : 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, &
475 : 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, &
476 : 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, &
477 : 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, &
478 : 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, 1.78650, &
479 : 1.78650, 1.78650, 1.78650, 1.78720, &
480 : 1.82961, 1.83258, 1.83149, 1.82748, 1.82224, &
481 : 1.81718, 1.81204, 1.80704, 1.80250, 1.79834, 1.79482, &
482 : 1.79214, 1.78843, 1.78601, 1.78434, 1.78322, 1.78248, &
483 : 1.78201, 1.78160, 1.78140, 1.78160, 1.78220, 1.78310, &
484 : 1.78380, 1.78390, 1.78400, 1.78400, 1.78400, 1.78400, &
485 : 1.78400, 1.78390, 1.78380, 1.78370, 1.78370, 1.78370, &
486 : 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, &
487 : 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, &
488 : 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, &
489 : 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, 1.78370, &
490 : 1.78370, 1.78400, 1.78450, &
491 : 1.82961, 1.83258, 1.83149, 1.82748, 1.82224, 1.81718, &
492 : 1.81204, 1.80704, 1.80250, 1.79834, 1.79482, 1.79214, &
493 : 1.78843, 1.78601, 1.78434, 1.78322, 1.78248, 1.78201, &
494 : 1.78150, 1.78070, 1.78010, 1.77890, 1.77790, 1.77730, &
495 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
496 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
497 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
498 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
499 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
500 : 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, 1.77720, &
501 : 1.77720, 1.77800/),shape=(/nwlt,4/))
502 :
503 : ! #####################################################################
504 : ! Defines wavelength dependent complex index of refraction for ice.
505 : ! Allowable wavelength range extends from 0.045 microns to 8.6 meter
506 : ! temperature dependence only considered beyond 167 microns.
507 : !
508 : ! interpolation is done n_r vs. log(xlam)
509 : ! n_r vs. t
510 : ! log(n_i) vs. log(xlam)
511 : ! log(n_i) vs. t
512 : !
513 : ! Stephen G. Warren - 1983
514 : ! Dept. of Atmospheric Sciences
515 : ! University of Washington
516 : ! Seattle, Wa 98195
517 : !
518 : ! Based on
519 : !
520 : ! Warren,S.G.,1984.
521 : ! Optical constants of ice from the ultraviolet to the microwave.
522 : ! Applied Optics,23,1206-1225
523 : !
524 : ! Reference temperatures are -1.0,-5.0,-20.0, and -60.0 deg C
525 : ! #####################################################################
526 :
527 13633725 : pi = acos(-1._wp)
528 13633725 : n_r = 0._wp
529 13633725 : n_i = 0._wp
530 13633725 : tk = t
531 :
532 : ! Convert frequency to wavelength (um)
533 13633725 : alam=3E5_wp/freq
534 13633725 : if((alam < wlmin) .or. (alam > wlmax)) then
535 0 : call errorMessage('FATAL ERROR(optics/optics_lib.f90:m_ice): wavelength out of bounds')
536 0 : return
537 : endif
538 :
539 13633725 : if (alam < cutice) then
540 : ! Region from 0.045 microns to 167.0 microns - no temperature depend
541 0 : do i=2,nwl
542 0 : if(alam < wl(i)) continue
543 : enddo
544 0 : x1 = log(wl(i-1))
545 0 : x2 = log(wl(i))
546 0 : y1 = tabre(i-1)
547 0 : y2 = tabre(i)
548 0 : x = log(alam)
549 0 : y = ((x-x1)*(y2-y1)/(x2-x1))+y1
550 0 : n_r = y
551 0 : y1 = log(abs(tabim(i-1)))
552 0 : y2 = log(abs(tabim(i)))
553 0 : y = ((x-x1)*(y2-y1)/(x2-x1))+y1
554 0 : n_i = exp(y)
555 : else
556 : ! Region from 167.0 microns to 8.6 meters - temperature dependence
557 13633725 : if(tk > temref(1)) tk=temref(1)
558 13633725 : if(tk < temref(4)) tk=temref(4)
559 30031194 : do i=2,4
560 30031194 : if(tk.ge.temref(i)) go to 12
561 : enddo
562 13633725 : 12 lt1 = i
563 13633725 : lt2 = i-1
564 299941950 : do i=2,nwlt
565 299941950 : if(alam.le.wlt(i)) go to 14
566 : enddo
567 13633725 : 14 x1 = log(wlt(i-1))
568 13633725 : x2 = log(wlt(i))
569 13633725 : y1 = tabret(i-1,lt1)
570 13633725 : y2 = tabret(i,lt1)
571 13633725 : x = log(alam)
572 13633725 : ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
573 13633725 : y1 = tabret(i-1,lt2)
574 13633725 : y2 = tabret(i,lt2)
575 13633725 : yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
576 13633725 : t1 = temref(lt1)
577 13633725 : t2 = temref(lt2)
578 13633725 : y = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
579 13633725 : n_r = y
580 13633725 : y1 = log(abs(tabimt(i-1,lt1)))
581 13633725 : y2 = log(abs(tabimt(i,lt1)))
582 13633725 : ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
583 13633725 : y1 = log(abs(tabimt(i-1,lt2)))
584 13633725 : y2 = log(abs(tabimt(i,lt2)))
585 13633725 : yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
586 13633725 : y = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
587 13633725 : n_i = exp(y)
588 : endif
589 : end subroutine m_ice
590 :
591 : ! ############################################################################
592 : ! subroutine MIEINT
593 : ! ############################################################################
594 2375789780 : Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error)
595 : ! ##########################################################################
596 : !
597 : ! General purpose Mie scattering routine for single particles
598 : ! Author: R Grainger 1990
599 : ! History:
600 : ! G Thomas, March 2005: Added calculation of Phase function and
601 : ! code to ensure correct calculation of backscatter coeficient
602 : ! Options/Extend_Source
603 : !
604 : ! ##########################################################################
605 : ! INPUTS
606 : integer, intent(in) :: &
607 : Inp
608 : real(wp),intent(in) :: &
609 : Dx !
610 : real(wp),intent(in),dimension(Inp) :: &
611 : Dqv
612 : Complex(wp),intent(in) :: &
613 : SCm!
614 :
615 : ! OUTPUTS
616 : Complex(wp),intent(out),dimension(InP) :: &
617 : Xs1, & !
618 : Xs2 !
619 : real(wp),intent(out) :: &
620 : Dqxt, & !
621 : Dqsc, & !
622 : Dg, & !
623 : Dbsc !
624 : real(wp),intent(out),dimension(InP) :: &
625 : DPh
626 : integer :: &
627 : Error !!
628 :
629 : ! PARAMETERS
630 : Integer,parameter :: &
631 : Imaxx = 12000, & !
632 : Itermax = 30000, & ! Must be large enough to cope with the
633 : ! largest possible nmx = x * abs(scm) + 15
634 : ! or nmx = Dx + 4.05*Dx**(1./3.) + 2.0
635 : Imaxnp = 10000 ! Change this as required
636 : Real(wp),parameter :: &
637 : RIMax=2.5, & ! Largest real part of refractive index
638 : IRIMax = -2 ! Largest imaginary part of refractive index
639 :
640 : ! Internal variables
641 : Integer :: I, NStop, NmX, N, Inp2
642 : Real(wp) :: Chi,Chi0,Chi1,APsi,APsi0,APsi1,Psi,Psi0,Psi1
643 : Real(wp),dimension(Imaxnp) :: Pi0,Pi1,Taun
644 : Complex(wp) :: Ir,Cm,A,ANM1,APB,B,BNM1,AMB,Xi,Xi0,Xi1,Y
645 : Complex(wp),dimension(Itermax) :: D
646 : Complex(wp),dimension(Imaxnp) :: Sp,Sm!
647 :
648 : ! ACCELERATOR VARIABLES
649 : Integer :: Tnp1,Tnm1
650 : Real(wp) :: Dn, Rnx,Turbo,A2
651 : real(wp),dimension(Imaxnp) :: S,T
652 : Complex(wp) :: A1
653 :
654 2375789780 : If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
655 0 : Error = 1
656 0 : Return
657 : EndIf
658 2375789780 : Cm = SCm
659 2375789780 : Ir = 1 / Cm
660 2375789780 : Y = Dx * Cm
661 2375789780 : If (Dx.Lt.0.02) Then
662 : NStop = 2
663 : Else
664 1204750670 : If (Dx.Le.8.0) Then
665 1176117184 : NStop = Dx + 4.00*Dx**(1./3.) + 2.0
666 : Else
667 28633486 : If (Dx.Lt. 4200.0) Then
668 28633486 : NStop = Dx + 4.05*Dx**(1./3.) + 2.0
669 : Else
670 0 : NStop = Dx + 4.00*Dx**(1./3.) + 2.0
671 : End If
672 : End If
673 : End If
674 2375789780 : NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
675 2375789780 : If (Nmx .gt. Itermax) then
676 0 : Error = 1
677 0 : Return
678 : End If
679 2375789780 : Inp2 = Inp+1
680 : !ds D(NmX) = cmplx(0,0,Kind=Kind(0d0))
681 2375789780 : D(NmX) = cmplx(0,0,Kind=wp)
682 46291719871 : Do N = Nmx-1,1,-1
683 43915930091 : A1 = (N+1) / Y
684 46291719871 : D(N) = A1 - 1/(A1+D(N+1))
685 : End Do
686 7127369340 : Do I =1,Inp2
687 4751579560 : Sm(I) = cmplx(0,0,Kind=wp)
688 : !ds Sm(I) = cmplx(0,0,Kind=Kind(0d0))
689 4751579560 : Sp(I) = cmplx(0,0,Kind=wp)
690 : !ds Sp(I) = cmplx(0,0,Kind=Kind(0d0))
691 4751579560 : Pi0(I) = 0
692 7127369340 : Pi1(I) = 1
693 : End Do
694 2375789780 : Psi0 = Cos(Dx)
695 2375789780 : Psi1 = Sin(Dx)
696 2375789780 : Chi0 =-Sin(Dx)
697 2375789780 : Chi1 = Cos(Dx)
698 2375789780 : APsi0 = Psi0
699 2375789780 : APsi1 = Psi1
700 2375789780 : Xi0 = cmplx(APsi0,Chi0,Kind=wp)
701 : !ds Xi0 = cmplx(APsi0,Chi0,Kind=Kind(0d0))
702 2375789780 : Xi1 = cmplx(APsi1,Chi1,Kind=wp)
703 : !ds Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
704 2375789780 : Dg = 0
705 2375789780 : Dqsc = 0
706 2375789780 : Dqxt = 0
707 2375789780 : Tnp1 = 1
708 12264930518 : Do N = 1,Nstop
709 9889140738 : DN = N
710 9889140738 : Tnp1 = Tnp1 + 2
711 9889140738 : Tnm1 = Tnp1 - 2
712 9889140738 : A2 = Tnp1 / (DN*(DN+1._wp))
713 : !ds A2 = Tnp1 / (DN*(DN+1D0))
714 9889140738 : Turbo = (DN+1._wp) / DN
715 : !ds Turbo = (DN+1D0) / DN
716 9889140738 : Rnx = DN/Dx
717 9889140738 : Psi = Tnm1*Psi1/Dx - Psi0
718 : !ds Psi = Dble(Tnm1)*Psi1/Dx - Psi0
719 9889140738 : APsi = Psi
720 9889140738 : Chi = Tnm1*Chi1/Dx - Chi0
721 9889140738 : Xi = cmplx(APsi,Chi,Kind=wp)
722 : !ds Xi = cmplx(APsi,Chi,Kind=Kind(0d0))
723 9889140738 : A = ((D(N)*Ir+Rnx)*APsi-APsi1) / ((D(N)*Ir+Rnx)* Xi- Xi1)
724 9889140738 : B = ((D(N)*Cm+Rnx)*APsi-APsi1) / ((D(N)*Cm+Rnx)* Xi- Xi1)
725 9889140738 : Dqxt = Tnp1*(A + B)+ Dqxt
726 : !ds Dqxt = Tnp1 * Dble(A + B) + Dqxt
727 9889140738 : Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
728 9889140738 : If (N.Gt.1) then
729 7513350958 : Dg = Dg + (dN*dN - 1) * (ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 *(ANM1*Conjg(BNM1)) / (dN*dN - dN)
730 : !ds Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
731 : End If
732 9889140738 : Anm1 = A
733 9889140738 : Bnm1 = B
734 9889140738 : APB = A2 * (A + B)
735 9889140738 : AMB = A2 * (A - B)
736 29667422214 : Do I = 1,Inp2
737 19778281476 : If (I.GT.Inp) Then
738 9889140738 : S(I) = -Pi1(I)
739 : Else
740 9889140738 : S(I) = Dqv(I) * Pi1(I)
741 : End If
742 19778281476 : T(I) = S(I) - Pi0(I)
743 19778281476 : Taun(I) = N*T(I) - Pi0(I)
744 19778281476 : Sp(I) = APB * (Pi1(I) + Taun(I)) + Sp(I)
745 19778281476 : Sm(I) = AMB * (Pi1(I) - Taun(I)) + Sm(I)
746 19778281476 : Pi0(I) = Pi1(I)
747 29667422214 : Pi1(I) = S(I) + T(I)*Turbo
748 : End Do
749 9889140738 : Psi0 = Psi1
750 9889140738 : Psi1 = Psi
751 9889140738 : Apsi1 = Psi1
752 9889140738 : Chi0 = Chi1
753 9889140738 : Chi1 = Chi
754 12264930518 : Xi1 = cmplx(APsi1,Chi1,Kind=wp)
755 : !ds Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
756 : End Do
757 :
758 2375789780 : If (Dg .GT.0) Dg = 2 * Dg / Dqsc
759 2375789780 : Dqsc = 2 * Dqsc / Dx**2
760 2375789780 : Dqxt = 2 * Dqxt / Dx**2
761 4751579560 : Do I = 1,Inp
762 2375789780 : Xs1(I) = (Sp(I)+Sm(I)) / 2
763 2375789780 : Xs2(I) = (Sp(I)-Sm(I)) / 2
764 4751579560 : Dph(I) = 2 * (Xs1(I)*Conjg(Xs1(I)) + Xs2(I)*Conjg(Xs2(I))) / (Dx**2 * Dqsc)
765 : !ds Dph(I) = 2 * Dble(Xs1(I)*Conjg(Xs1(I)) + Xs2(I)*Conjg(Xs2(I))) / (Dx**2 * Dqsc)
766 : End Do
767 2375789780 : Dbsc = 4 * Abs(( (Sp(Inp2)+Sm(Inp2))/2 )**2) / Dx**2
768 2375789780 : Error = 0
769 2375789780 : Return
770 : End subroutine MieInt
771 : end module optics_lib
|
