source: trunk/htdocs/alexei_aercomp.html @ 1

Last change on this file since 1 was 1, checked in by cbipsl, 18 years ago

Geisa inital import
File size: 13.3 KB
Line 
1<!doctype html public "-//w3c//dtd html 4.0 transitional//en">
2<html>
3<head>
4   <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
5   <meta name="Author" content="Alexey Chursin">
6   <meta name="GENERATOR" content="Mozilla/4.7 [en] (Win95; I) [Netscape]">
7   <title>The Aecomp.f program source</title>
8</head>
9<body>
10!*****&nbsp;&nbsp;&nbsp; PROGRAM&nbsp; FOR COMPUTATION
11OF THE INTEGRATED OPTICAL CHARACTERISTICS ***
12<br>!*****&nbsp;&nbsp;&nbsp; OF ATMOSPHERIC AEROSOLS&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
13*****
14<br>!***&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ( NC_MAX&nbsp;&nbsp;
15- MAX. NUMBER OF COMPONENTS&nbsp; )&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
16**
17<br>!***&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ( NP&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
18- MAX. NUMBER OF ANGULAR POINTS&nbsp; )&nbsp;&nbsp;&nbsp; **
19<br>!***&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ( NDV_MAX&nbsp;
20- MAX. NUMBER OF WAVELENGTHS&nbsp; )&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
21**
22<br>!*****************************************************************
23<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; PARAMETER (NC_MAX=5,NP=204,NDV_MAX=70,
24<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
25pi=3.14159265358D0)
26<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; REAL*4&nbsp;&nbsp;&nbsp; dlm,lm1,lm3,ex1,ex3,w1,w3,
27<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
28ex(NC_MAX),w(NC_MAX),smu(NC_MAX),s(NC_MAX),
29<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
30ex0,w0,smu0,s0
31<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; REAL*8&nbsp;&nbsp;&nbsp; mu,ug,ug_n,mu_n,dmu
32<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; INTEGER&nbsp;&nbsp; n_ug
33<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CHARACTER s_dat*12,ind_dat*12,
34<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
35name_m0*12,name_ug0*12,
36<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
37name_m*12, name_ug*12,
38<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
39nc_ind*12,buf*80
40<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DIMENSION n_ug(NC_MAX),
41<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
42d1(NP),d2(NC_MAX,NP),d3(NP),ug(NP),mu(NP),dv(NDV_MAX),
43<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
44d0(NP),si(NC_MAX),n_m(NC_MAX),
45<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
46ug_n(NC_MAX,NP),mu_n(NC_MAX,NP),dv_n(NC_MAX,NDV_MAX),
47<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
48name_m(NC_MAX), name_ug(NC_MAX),
49<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
50nc_ind(NC_MAX)
51<br>c *****
52<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN(11,file = 'inp_comp.dat')
53<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ (11,43) n_comp,n_m0,name_m0,n_ug0,name_ug0,
54<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
55s_dat,ind_dat,buf
56<br>43&nbsp;&nbsp;&nbsp; FORMAT(I6/I6/A12/I6/A12/A12/A12/A50)
57<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; PRINT 50, n_comp,n_m0,name_m0,n_ug0,name_ug0,s_dat,ind_dat
58<br>50&nbsp;&nbsp;&nbsp; FORMAT (/
59<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
60NUMBER OF AEROSOL COMPONENTS: ',I3/
61<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
62NUMBER OF WAVELENGTHS: ',I3/
63<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
64NAME&nbsp; OF THE FILE WITH WAVELENGTHS: ',A/
65<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; TOTAL
66NUMBER OF ANGULAR NODAL POINTS: ',I3/
67<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' ** NAME&nbsp; OF THE FILE OF ANGULAR NODAL
68POINTS: ',A/
69<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp; NAME&nbsp; OF THE OUTPUT
70FILE FOR EXTINCTION: ',A/
71<br>&nbsp;&nbsp;&nbsp;&nbsp; + ' **&nbsp;&nbsp; NAME&nbsp; OF THE OUTPUT
72FILE FOR PHASE FUNCTION: ',A/)
73<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; pause
74<br>c *************************************************
75<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN&nbsp; (41,file=s_dat)
76<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN&nbsp; (21,file=name_m0)
77<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ(21,*)(dv(i),i=1,n_m0)
78<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CLOSE (21)
79<br>c *************************************************
80<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN&nbsp; (21,file=name_ug0)
81<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ&nbsp; (21,*) (ug(i),i=1,n_ug0)
82<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CLOSE (21)
83<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i1=1,n_ug0
84<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; mu(i1)=dcos(pi/180.d0*ug(i1))
85<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
86<br>c *************************************************
87<br>c ******&nbsp;&nbsp;&nbsp; BEGINNING OF COMPONENT LOOP&nbsp;&nbsp;&nbsp;
88*************
89<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i0 = 1,n_comp
90<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ (11,111) n_m(i0),nc_ind(i0),si(i0),
91<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
92name_m(i0),n_ug(i0),name_ug(i0),buf
93<br>111&nbsp;&nbsp; FORMAT(I6/A12/E14.6/A12/I6/A12/A50)
94<br>&nbsp; Write(*,102)
95<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; *****&nbsp; COMPONENT NUMBER:&nbsp;&nbsp;
96', I0, '&nbsp; *****',
97<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp; NAME&nbsp;
98OF THE INPUT FILE FOR PHASE FUNCTION: ',
99<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; nc_ind(i0),
100<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
101NUMBER OF PARTICLES OF THIS KIND: ',si(i0),
102<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
103NUMBER OF WAVELENGTHS: ',n_m(i0),
104<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
105NAME OF THE FILE FOR WAVELENGTHS: ',name_m(i0),
106<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp;&nbsp;&nbsp;
107TOTAL NUMBER OF ANGULAR NODAL POINTS: ',n_ug(i0),
108<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp; '&nbsp; **&nbsp;&nbsp;&nbsp;&nbsp;
109NAME&nbsp; OF THE FILE FOR ANGULAR POINTS: ',name_ug(i0)
110<br>102&nbsp;&nbsp;&nbsp;&nbsp; FORMAT(A,I3,A/A48,A,/A48,E12.4/A48,I6/A48,A12/A48,i6/A48,A12)
111<br>c ******************************************************************
112<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; pause
113<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN&nbsp; (21,file=name_m(i0))
114<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ(21,*)(dv_n(i0,i),i=1,n_m(i0))
115<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CLOSE (21)
116<br>c ***********
117<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN&nbsp; (21,file=name_ug(i0))
118<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ&nbsp; (21,*) (ug_n(i0,i),i=1,n_ug(i0))
119<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CLOSE (21)
120<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i1=1,n_ug(i0)
121<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; mu_n(i0,i1)=dcos(pi/180.d0*ug_n(i0,i1))
122<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
123<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
124<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; PRINT *, ' ***************************************************'
125<br>c ***************************************************************
126<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN (66,FILE = ind_dat, ACCESS='DIRECT',
127<br>&nbsp;&nbsp;&nbsp;&nbsp; +FORM = 'UNFORMATTED',RECL=4*n_ug0+24 ) !&nbsp;
128NN+lm+1+...+n_ug
129<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO il=1,n_m0
130<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dlm=dv(il)
131<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ex0=0.
132<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; e_sct =0.
133<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i=1,n_ug0
134<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; d0(i)=0.
135<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
136<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i0=1,n_comp
137<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; iv =1
138<br>20&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; iv=iv+1
139<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; IF (iv .GT.
140n_m(i0) .OR. dlm .LT. dv_n(i0,1))THEN
141<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
142PRINT *, '** ERROR OF LENGTH-WAVE IN ',name_m(i0)
143<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
144STOP
145<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END IF
146<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; IF (dlm .GT.
147dv_n(i0,iv)) GOTO 20
148<br>c ********
149<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OPEN (61,FILE
150= nc_ind(i0), ACCESS='DIRECT',
151<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp; FORM = 'UNFORMATTED',RECL=4*n_ug(i0)+24
152) !&nbsp; NN+lm+1+...+n_ug
153<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ(61,REC=iv-1)q1,lm1,ex1,w1,smu1,s1,(d1(I),i=1,n_ug(i0))
154<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; READ(61,REC=iv)&nbsp;
155q3,lm3,ex3,w3,smu3,s3,(d3(I),i=1,n_ug(i0))
156<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; CLOSE(61)
157<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i1=1,n_ug(i0)
158<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
159d2(i0,i1)=d1(i1)+(d3(i1)-d1(i1))/(lm3-lm1)*(dlm-lm1)
160<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
161<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ex(i0)=ex1+(ex3-ex1)/(lm3-lm1)*(dlm-lm1)
162<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; w(i0) = w1+(
163w3- w1)/(lm3-lm1)*(dlm-lm1)
164<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; smu(i0) = smu1+(smu3-
165smu1)/(lm3-lm1)*(dlm-lm1)
166<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; s(i0) = s1+(
167s3- s1)/(lm3-lm1)*(dlm-lm1)
168<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ex0=ex0+ex(i0)*si(i0)
169<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; e_sct=e_sct+ex(i0)*w(i0)*si(i0)
170<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
171<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; w0= e_sct/ex0
172<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i1=1,n_ug0
173<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; d0(i1)=0.
174<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO i0=1,n_comp
175<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
176iv =1
177<br>25&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; iv=iv+1
178<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
179IF(mu(i1).GE.mu_n(i0,iv).AND.mu(i1).LE.mu_n(i0,iv-1))THEN
180<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
181ELSE
182<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
183GOTO 25
184<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
185END IF
186<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
187dmu=mu(i1) - mu_n(i0,iv-1)
188<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
189dd2=d2(i0,iv)-d2(i0,iv-1)
190<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
191d10=d2(i0,iv-1)+dd2/
192<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
193(mu_n(i0,iv)-mu_n(i0,iv-1))* dmu
194<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
195dd2=d10*ex(i0)*w(i0)*si(i0)/e_sct
196<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
197d0(i1)=d0(i1)+ dd2
198<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
199<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
200<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; s0=0
201<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; smu0=0
202<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO j=1,n_ug0-1
203<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; s0&nbsp;
204=s0&nbsp; +(d0(j)+d0(j+1))/2.*(mu(j)-mu(j+1))
205<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; smu0=smu0+d0(j)*mu(j)*(mu(j)-mu(j+1))+(d0(j+1)-d0(j))/2.*
206<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
207(mu(j)-mu(j+1))
208<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
209<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; q1=il
210<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; smu0=smu0/s0
211<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; s0=2.*pi*s0
212<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; WRITE(66,REC=il)q1,dlm,ex0,w0,smu0,s0,(d0(I),i=1,n_ug0)
213<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; WRITE (41,'(I3,2x,F6.3,2x,3(E12.5,2x),F6.3)')
214<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
215il,dlm,ex0,w0,smu0,s0
216<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; WRITE (* ,'(I3,2x,F6.3,2x,3(E12.5,2x),F6.3)')
217<br>&nbsp;&nbsp;&nbsp;&nbsp; +&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
218il,dlm,ex0,w0,smu0,s0
219<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END DO
220<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END
221</body>
222</html>
Note: See TracBrowser for help on using the repository browser.