source: LMDZ4/trunk/libf/cosp/optics_lib.F90 @ 5295

Last change on this file since 5295 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
File size: 35.0 KB
RevLine 
[1262]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
Note: See TracBrowser for help on using the repository browser.