source: LMDZ5/trunk/libf/phylmd/cosp/optics_lib.F90 @ 3740

Last change on this file since 3740 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
RevLine 
[2428]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 $
[1262]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! ----------------------------------------------------------------------------
[2428]19  subroutine m_wat(freq, tk, n_r, n_i)
[1262]20  implicit none
21! 
22! Purpose:
23!   compute complex index of refraction of liquid water
24!
25! Inputs:
26!   [freq]    frequency (GHz)
[2428]27!   [tk]       temperature (K)
[1262]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 -----
[2428]40  real*8, intent(in) :: freq,tk
[1262]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
[2428]49  real*8 tc
[1262]50  complex*16 e_comp, sq
51
[2428]52  tc = tk - 273.15
53
[1262]54  ld = 100.*2.99792458E8/(freq*1E9)
[2428]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.021647*tc-0.00131198*tc**2
58  a = -(16.8129/(tc+273.))+0.0609265
59  ls = 0.00033836*exp(2513.98/(tc+273.))
[1262]60  sg = 12.5664E8
61
62  tm1 = (ls/ld)**(1-a)
63  pi = acos(-1.D0)
64  cos1 = cos(0.5*a*pi)
65  sin1 = sin(0.5*a*pi)
66
67  e_r = ei + (((es-ei)*(1.+tm1*sin1))/(1.+2*tm1*sin1+tm1**2))
68  e_i = (((es-ei)*tm1*cos1)/(1.+2*tm1*sin1+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)
[2428]91!   [t]       temperature (K)
[1262]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
[2428]113  real*8 :: alam,cutice,pi,t1,t2,wlmax,wlmin, &
114            x,x1,x2,y,y1,y2,ylo,yhi,tk
[1262]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
[2428]509  tk = t
510
[1262]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.05*Dx**(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.00*Dx**(1./3.) + 2.0
659         Else
660            If (Dx.Lt. 4200.0) Then
661               NStop = Dx + 4.05*Dx**(1./3.) + 2.0
662            Else
663               NStop = Dx + 4.00*Dx**(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 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
711         If (N.Gt.1) then
[2428]712         Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
[1262]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) / Dx**2
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.