source: LMDZ6/trunk/libf/phylmdiso/cosp2/optics_lib.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 40.2 KB
Line 
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! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34module optics_lib
35  USE COSP_KINDS,     ONLY: wp
36  use mod_cosp_error, ONLY: errorMessage
37  implicit none
38
39contains
40
41  ! ##############################################################################
42  !                           subroutine M_WAT
43  ! ##############################################################################
44  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    tc = tk - 273.15_wp
80
81    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         -2.8E-8_wp*(tc-25._wp)**3))
84    ei = 5.27137_wp+0.021647_wp*tc-0.00131198_wp*tc**2
85    a = -(16.8129_wp/(tc+273._wp))+0.0609265_wp
86    ls = 0.00033836_wp*exp(2513.98_wp/(tc+273._wp))
87    sg = 12.5664E8_wp
88   
89    tm1 = (ls/ld)**(1-a)
90    pi = acos(-1._wp)
91    cos1 = cos(0.5_wp*a*pi)
92    sin1 = sin(0.5_wp*a*pi)
93   
94    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         +((sg*ld)/1.885E11_wp)
97   
98!ds    e_comp = cmplx(e_r,e_i,Kind=Kind(0d0))
99    e_comp = cmplx(e_r,e_i,Kind=wp)
100    sq = sqrt(e_comp)
101   
102    n_r = real(sq)
103    n_i = aimag(sq)     
104   
105    return
106  end subroutine m_wat
107
108  ! ############################################################################
109  !                           subroutine M_ICE
110  ! ############################################################################
111  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    pi  = acos(-1._wp)
528    n_r = 0._wp
529    n_i = 0._wp
530    tk  = t
531   
532    ! Convert frequency to wavelength (um)
533    alam=3E5_wp/freq
534    if((alam < wlmin) .or. (alam > wlmax)) then
535       call errorMessage('FATAL ERROR(optics/optics_lib.f90:m_ice): wavelength out of bounds')
536       return
537    endif
538   
539    if (alam < cutice) then
540       ! Region from 0.045 microns to 167.0 microns - no temperature depend
541       do i=2,nwl
542          if(alam < wl(i)) continue
543       enddo
544       x1  = log(wl(i-1))
545       x2  = log(wl(i))
546       y1  = tabre(i-1)
547       y2  = tabre(i)
548       x   = log(alam)
549       y   = ((x-x1)*(y2-y1)/(x2-x1))+y1
550       n_r = y
551       y1  = log(abs(tabim(i-1)))
552       y2  = log(abs(tabim(i)))
553       y   = ((x-x1)*(y2-y1)/(x2-x1))+y1
554       n_i = exp(y)   
555    else
556       ! Region from 167.0 microns to 8.6 meters - temperature dependence
557       if(tk > temref(1)) tk=temref(1)
558       if(tk < temref(4)) tk=temref(4)
559       do i=2,4
560          if(tk.ge.temref(i)) go to 12
561       enddo
56212     lt1 = i
563       lt2 = i-1
564       do i=2,nwlt
565          if(alam.le.wlt(i)) go to 14
566       enddo
56714     x1  = log(wlt(i-1))
568       x2  = log(wlt(i))
569       y1  = tabret(i-1,lt1)
570       y2  = tabret(i,lt1)
571       x   = log(alam)
572       ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
573       y1  = tabret(i-1,lt2)
574       y2  = tabret(i,lt2)
575       yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
576       t1  = temref(lt1)
577       t2  = temref(lt2)
578       y   = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
579       n_r = y
580       y1  = log(abs(tabimt(i-1,lt1)))
581       y2  = log(abs(tabimt(i,lt1)))
582       ylo = ((x-x1)*(y2-y1)/(x2-x1))+y1
583       y1  = log(abs(tabimt(i-1,lt2)))
584       y2  = log(abs(tabimt(i,lt2)))
585       yhi = ((x-x1)*(y2-y1)/(x2-x1))+y1
586       y   = ((tk-t1)*(yhi-ylo)/(t2-t1))+ylo
587       n_i = exp(y)
588    endif
589  end subroutine m_ice
590
591  ! ############################################################################
592  ! subroutine MIEINT
593  ! ############################################################################
594  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    If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
655       Error = 1
656       Return
657    EndIf
658    Cm = SCm
659    Ir = 1 / Cm
660    Y =  Dx * Cm
661    If (Dx.Lt.0.02) Then
662       NStop = 2
663    Else
664       If (Dx.Le.8.0) Then
665          NStop = Dx + 4.00*Dx**(1./3.) + 2.0
666       Else
667          If (Dx.Lt. 4200.0) Then
668             NStop = Dx + 4.05*Dx**(1./3.) + 2.0
669          Else
670             NStop = Dx + 4.00*Dx**(1./3.) + 2.0
671          End If
672       End If
673    End If
674    NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
675    If (Nmx .gt. Itermax) then
676       Error = 1
677       Return
678    End If
679    Inp2 = Inp+1
680!ds    D(NmX) = cmplx(0,0,Kind=Kind(0d0))
681    D(NmX) = cmplx(0,0,Kind=wp)
682    Do N = Nmx-1,1,-1
683       A1 = (N+1) / Y
684       D(N) = A1 - 1/(A1+D(N+1))
685    End Do
686    Do I =1,Inp2
687       Sm(I) = cmplx(0,0,Kind=wp)
688!ds       Sm(I) = cmplx(0,0,Kind=Kind(0d0))
689       Sp(I) = cmplx(0,0,Kind=wp)
690!ds       Sp(I) = cmplx(0,0,Kind=Kind(0d0))
691       Pi0(I) = 0
692       Pi1(I) = 1
693    End Do
694    Psi0 = Cos(Dx)
695    Psi1 = Sin(Dx)
696    Chi0 =-Sin(Dx)
697    Chi1 = Cos(Dx)
698    APsi0 = Psi0
699    APsi1 = Psi1
700    Xi0 = cmplx(APsi0,Chi0,Kind=wp)
701!ds    Xi0 = cmplx(APsi0,Chi0,Kind=Kind(0d0))
702    Xi1 = cmplx(APsi1,Chi1,Kind=wp)
703!ds    Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
704    Dg = 0
705    Dqsc = 0
706    Dqxt = 0
707    Tnp1 = 1
708    Do N = 1,Nstop
709       DN = N
710       Tnp1 = Tnp1 + 2
711       Tnm1 = Tnp1 - 2
712       A2 = Tnp1 / (DN*(DN+1._wp))
713!ds       A2 = Tnp1 / (DN*(DN+1D0))
714       Turbo = (DN+1._wp) / DN
715!ds       Turbo = (DN+1D0) / DN
716       Rnx = DN/Dx
717       Psi = Tnm1*Psi1/Dx - Psi0
718!ds       Psi = Dble(Tnm1)*Psi1/Dx - Psi0
719       APsi = Psi
720       Chi = Tnm1*Chi1/Dx       - Chi0
721       Xi = cmplx(APsi,Chi,Kind=wp)
722!ds       Xi = cmplx(APsi,Chi,Kind=Kind(0d0))
723       A = ((D(N)*Ir+Rnx)*APsi-APsi1) / ((D(N)*Ir+Rnx)*  Xi-  Xi1)
724       B = ((D(N)*Cm+Rnx)*APsi-APsi1) / ((D(N)*Cm+Rnx)*  Xi-  Xi1)
725       Dqxt = Tnp1*(A + B)+ Dqxt
726!ds       Dqxt = Tnp1 *      Dble(A + B)          + Dqxt
727       Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
728       If (N.Gt.1) then
729          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       Anm1 = A
733       Bnm1 = B
734       APB = A2 * (A + B)
735       AMB = A2 * (A - B)
736       Do I = 1,Inp2
737          If (I.GT.Inp) Then
738             S(I) = -Pi1(I)
739          Else
740             S(I) = Dqv(I) * Pi1(I)
741          End If
742          T(I) = S(I) - Pi0(I)
743          Taun(I) = N*T(I) - Pi0(I)
744          Sp(I) = APB * (Pi1(I) + Taun(I)) + Sp(I)
745          Sm(I) = AMB * (Pi1(I) - Taun(I)) + Sm(I)
746          Pi0(I) = Pi1(I)
747          Pi1(I) = S(I) + T(I)*Turbo
748       End Do
749       Psi0 = Psi1
750       Psi1 = Psi
751       Apsi1 = Psi1
752       Chi0 = Chi1
753       Chi1 = Chi
754       Xi1 = cmplx(APsi1,Chi1,Kind=wp)
755!ds       Xi1 = cmplx(APsi1,Chi1,Kind=Kind(0d0))
756    End Do
757
758    If (Dg .GT.0) Dg = 2 * Dg / Dqsc
759    Dqsc =  2 * Dqsc / Dx**2
760    Dqxt =  2 * Dqxt / Dx**2
761    Do I = 1,Inp
762       Xs1(I) = (Sp(I)+Sm(I)) / 2
763       Xs2(I) = (Sp(I)-Sm(I)) / 2
764       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    Dbsc = 4 * Abs(( (Sp(Inp2)+Sm(Inp2))/2 )**2) / Dx**2
768    Error = 0
769    Return
770  End subroutine MieInt
771end module optics_lib
Note: See TracBrowser for help on using the repository browser.