CDS-ESPRI: aercomp

!*****    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