!***** PROGRAM FOR COMPUTATION
OF THE INTEGRATED OPTICAL CHARACTERISTICS ***
!***** OF ATMOSPHERIC AEROSOLS
*****
!*** ( NC_MAX
- MAX. NUMBER OF COMPONENTS )
**
!*** ( NP
- MAX. NUMBER OF ANGULAR POINTS ) **
!*** ( NDV_MAX
- MAX. NUMBER OF WAVELENGTHS )
**
!*****************************************************************
PARAMETER (NC_MAX=5,NP=204,NDV_MAX=70,
+
pi=3.14159265358D0)
REAL*4 dlm,lm1,lm3,ex1,ex3,w1,w3,
+
ex(NC_MAX),w(NC_MAX),smu(NC_MAX),s(NC_MAX),
+
ex0,w0,smu0,s0
REAL*8 mu,ug,ug_n,mu_n,dmu
INTEGER n_ug
CHARACTER s_dat*12,ind_dat*12,
+
name_m0*12,name_ug0*12,
+
name_m*12, name_ug*12,
+
nc_ind*12,buf*80
DIMENSION n_ug(NC_MAX),
+
d1(NP),d2(NC_MAX,NP),d3(NP),ug(NP),mu(NP),dv(NDV_MAX),
+
d0(NP),si(NC_MAX),n_m(NC_MAX),
+
ug_n(NC_MAX,NP),mu_n(NC_MAX,NP),dv_n(NC_MAX,NDV_MAX),
+
name_m(NC_MAX), name_ug(NC_MAX),
+
nc_ind(NC_MAX)
c *****
OPEN(11,file = 'inp_comp.dat')
READ (11,43) n_comp,n_m0,name_m0,n_ug0,name_ug0,
+
s_dat,ind_dat,buf
43 FORMAT(I6/I6/A12/I6/A12/A12/A12/A50)
PRINT 50, n_comp,n_m0,name_m0,n_ug0,name_ug0,s_dat,ind_dat
50 FORMAT (/
+ ' **
NUMBER OF AEROSOL COMPONENTS: ',I3/
+ ' **
NUMBER OF WAVELENGTHS: ',I3/
+ ' **
NAME OF THE FILE WITH WAVELENGTHS: ',A/
+ ' ** TOTAL
NUMBER OF ANGULAR NODAL POINTS: ',I3/
+ ' ** NAME OF THE FILE OF ANGULAR NODAL
POINTS: ',A/
+ ' ** NAME OF THE OUTPUT
FILE FOR EXTINCTION: ',A/
+ ' ** NAME OF THE OUTPUT
FILE FOR PHASE FUNCTION: ',A/)
pause
c *************************************************
OPEN (41,file=s_dat)
OPEN (21,file=name_m0)
READ(21,*)(dv(i),i=1,n_m0)
CLOSE (21)
c *************************************************
OPEN (21,file=name_ug0)
READ (21,*) (ug(i),i=1,n_ug0)
CLOSE (21)
DO i1=1,n_ug0
mu(i1)=dcos(pi/180.d0*ug(i1))
END DO
c *************************************************
c ****** BEGINNING OF COMPONENT LOOP
*************
DO i0 = 1,n_comp
READ (11,111) n_m(i0),nc_ind(i0),si(i0),
+
name_m(i0),n_ug(i0),name_ug(i0),buf
111 FORMAT(I6/A12/E14.6/A12/I6/A12/A50)
Write(*,102)
+ ' ***** COMPONENT NUMBER:
', I0, ' *****',
+ ' ** NAME
OF THE INPUT FILE FOR PHASE FUNCTION: ',
+ nc_ind(i0),
+ ' **
NUMBER OF PARTICLES OF THIS KIND: ',si(i0),
+ ' **
NUMBER OF WAVELENGTHS: ',n_m(i0),
+ ' **
NAME OF THE FILE FOR WAVELENGTHS: ',name_m(i0),
+ ' **
TOTAL NUMBER OF ANGULAR NODAL POINTS: ',n_ug(i0),
+ ' **
NAME OF THE FILE FOR ANGULAR POINTS: ',name_ug(i0)
102 FORMAT(A,I3,A/A48,A,/A48,E12.4/A48,I6/A48,A12/A48,i6/A48,A12)
c ******************************************************************
pause
OPEN (21,file=name_m(i0))
READ(21,*)(dv_n(i0,i),i=1,n_m(i0))
CLOSE (21)
c ***********
OPEN (21,file=name_ug(i0))
READ (21,*) (ug_n(i0,i),i=1,n_ug(i0))
CLOSE (21)
DO i1=1,n_ug(i0)
mu_n(i0,i1)=dcos(pi/180.d0*ug_n(i0,i1))
END DO
END DO
PRINT *, ' ***************************************************'
c ***************************************************************
OPEN (66,FILE = ind_dat, ACCESS='DIRECT',
+FORM = 'UNFORMATTED',RECL=4*n_ug0+24 ) !
NN+lm+1+...+n_ug
DO il=1,n_m0
dlm=dv(il)
ex0=0.
e_sct =0.
DO i=1,n_ug0
d0(i)=0.
END DO
DO i0=1,n_comp
iv =1
20 iv=iv+1
IF (iv .GT.
n_m(i0) .OR. dlm .LT. dv_n(i0,1))THEN
PRINT *, '** ERROR OF LENGTH-WAVE IN ',name_m(i0)
STOP
END IF
IF (dlm .GT.
dv_n(i0,iv)) GOTO 20
c ********
OPEN (61,FILE
= nc_ind(i0), ACCESS='DIRECT',
+ FORM = 'UNFORMATTED',RECL=4*n_ug(i0)+24
) ! NN+lm+1+...+n_ug
READ(61,REC=iv-1)q1,lm1,ex1,w1,smu1,s1,(d1(I),i=1,n_ug(i0))
READ(61,REC=iv)
q3,lm3,ex3,w3,smu3,s3,(d3(I),i=1,n_ug(i0))
CLOSE(61)
DO i1=1,n_ug(i0)
d2(i0,i1)=d1(i1)+(d3(i1)-d1(i1))/(lm3-lm1)*(dlm-lm1)
END DO
ex(i0)=ex1+(ex3-ex1)/(lm3-lm1)*(dlm-lm1)
w(i0) = w1+(
w3- w1)/(lm3-lm1)*(dlm-lm1)
smu(i0) = smu1+(smu3-
smu1)/(lm3-lm1)*(dlm-lm1)
s(i0) = s1+(
s3- s1)/(lm3-lm1)*(dlm-lm1)
ex0=ex0+ex(i0)*si(i0)
e_sct=e_sct+ex(i0)*w(i0)*si(i0)
END DO
w0= e_sct/ex0
DO i1=1,n_ug0
d0(i1)=0.
DO i0=1,n_comp
iv =1
25 iv=iv+1
IF(mu(i1).GE.mu_n(i0,iv).AND.mu(i1).LE.mu_n(i0,iv-1))THEN
ELSE
GOTO 25
END IF
dmu=mu(i1) - mu_n(i0,iv-1)
dd2=d2(i0,iv)-d2(i0,iv-1)
d10=d2(i0,iv-1)+dd2/
+
(mu_n(i0,iv)-mu_n(i0,iv-1))* dmu
dd2=d10*ex(i0)*w(i0)*si(i0)/e_sct
d0(i1)=d0(i1)+ dd2
END DO
END DO
s0=0
smu0=0
DO j=1,n_ug0-1
s0
=s0 +(d0(j)+d0(j+1))/2.*(mu(j)-mu(j+1))
smu0=smu0+d0(j)*mu(j)*(mu(j)-mu(j+1))+(d0(j+1)-d0(j))/2.*
+
(mu(j)-mu(j+1))
END DO
q1=il
smu0=smu0/s0
s0=2.*pi*s0
WRITE(66,REC=il)q1,dlm,ex0,w0,smu0,s0,(d0(I),i=1,n_ug0)
WRITE (41,'(I3,2x,F6.3,2x,3(E12.5,2x),F6.3)')
+
il,dlm,ex0,w0,smu0,s0
WRITE (* ,'(I3,2x,F6.3,2x,3(E12.5,2x),F6.3)')
+
il,dlm,ex0,w0,smu0,s0
END DO
END