Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

optics_lib.F90

Last change on this file was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 35.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format