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