Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  Main
Language: FORTRAN
Code:
C   THIS CODE COMPUTES THE BIDIRECTIONAL REFLECTION FUNCTION
C   FOR A SEMI-INFINITE PARTICULATE SLAB USING THE ARRAY
C   OF PRECALCULATED FOURIER COMPONENTS OF THE REFLECTION
C   FUNCTION SPECIFIED AT THE DIVISION NODES OF A QUADRATURE
C   FORMULA ON THE INTERVAL (0, 1) (FILE 11).

C   INPUT PARAMETERS:
C
C   MU0 - COSINE OF THE ILLUMINATION ZENITH ANGLE 
C   MU - NMU VALUES OF THE COSINE OF THE REFLECTION ZENITH
C        ANGLE EVENLY DISTRIBUTED OVER THE INTERVAL (0, 1]
C        (THE SMALLEST MU VALUE IS 0.01)
C   PHI - NPHI VALUES OF THE AZIMUTH ANGLE EVENLY DISTRIBUTED
C   OVER THE INTERVAL [0, 180] DEGREES

C   OUTPUT PARAMETERS:
C
C   RO - ARRAY OF THE RESPECTIVE BRF VALUES (FILE 14)

C   MESSAGE 'NGS.NE.NG, EXECUTION TERMINATED' MEANS THAT
C   THE NUMBER OF QUADRATURE POINTS IN THE INTERPOLATION
C   PROGRAM IS NOT EQUAL TO THE RESPECTIVE NUMBER IN THE
C   PROGRAM COMPUTING THE FOURIER COMPONENTS OF THE 
C   REFLECTION FUNCTION AND MUST BE CHANGED ACCORDINGLY.
  
C   FOR MORE INFORMATION SEE THE DESCRIPTION OF THE CODE COMPUTING
C   THE FOURIER COMPONENTS OF THE REFLECTION FUNCTION.

C   WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
C   IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
C   INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
C   VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
C   THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
C   ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM. 


      PARAMETER (LMAX1=700, NMU=41, NPHI=91, NG=100)
      IMPLICIT REAL*8 (A-H,O-Z)                        
      REAL*8 X(NG),MU(NMU),PHI(NPHI),                  
     *       AL1(LMAX1),MU0,RI(NG),
     *       BB(NG),CC(NG),DD(NG),Y(NG),
     *       R(NG,NG),RR(NMU,LMAX1),RO(NMU,NPHI)
      COMMON /COEFF/ AL1,MMAX1
                                                                        
      open (11,file='refl.write')
      open (14,file='interp.write')

      MU0=0.7D0
      DO I=1,NMU
         MU(I)=DFLOAT(I-1)/DFLOAT(NMU-1)
      ENDDO
      IF (MU(1).LT.0.01D0) MU(1)=0.01D0
      DO I=1,NPHI
         PHI(I)=180D0*DFLOAT(I-1)/DFLOAT(NPHI-1)
      ENDDO
      READ (11,*) NSEP

C***  READ IN QUADRATURE DIVISION POINTS
      READ (11,1501) NGS
      IF (NG.NE.NGS) WRITE (6,1001)
      IF (NG.NE.NGS) STOP      
 1001 FORMAT ('NGS.NE.NG, EXECUTION TERMINATED')
      READ (11,1502) X
 1501 FORMAT (I5,5D15.8)
 1502 FORMAT (6D15.8)
 1503 FORMAT (3I5)

C***  READ IN SINGLE-SCATTERING ALBEDO AND LEGENDRE
C***  EXPANSION COEFFICIENTS
      READ (11,1501) MMAX1,ALB
      IF (MMAX1.GT.LMAX1) WRITE (6,1002)
 1002 FORMAT ('INCREASE LMAX1 IN ALL PARAMETER STATEMENTS!')
      IF (MMAX1.GT.LMAX1) STOP
      READ (11,1502) (AL1(M1), M1=1,MMAX1)

C***  WRITE MU0, MU, AND PHI VALUES  
      WRITE (14,1502) MU0
      WRITE (14,1503) NMU, NPHI
      WRITE (14,1502) MU
      WRITE (14,1502) PHI
      DO I=1,NMU                                                   
         IF (DABS(MU(I)-1D0).LT.1D-9) MU(I)=0.999999999D0
      ENDDO                                                            

      DO M1=1,MMAX1                                                 
         READ (11,*) IFLAG

C***  READ IN (M1-1)th FOURIER COMPONENT OF THE REFLECTION FUNCTION
         DO I=1,NG
            READ (11,1502) (R(I,J), J=1,I)
            DO J=1,I
               R(J,I)=R(I,J)
            ENDDO
         ENDDO
         MMAX2=M1

         DO I=1,NG                                                
            DO J=1,NG                                             
               Y(J)=R(I,J)*X(J)*X(I)            
            ENDDO                              
            CALL SPLINE (NG,X,Y,BB,CC,DD)     
            RI(I)=SEVAL(NG,MU0,X,Y,BB,CC,DD)     
         ENDDO
         CALL SPLINE (NG,X,RI,BB,CC,DD)
         DO K=1,NMU
            RR(K,M1)=SEVAL(NG,MU(K),X,RI,BB,CC,DD)/MU(K)
         ENDDO
         IF (IFLAG.EQ.0) GO TO 100
      ENDDO                       
  100 CONTINUE
      P=DACOS(-1D0)/180D0
      DO I=1,NMU
         DO J=1,NPHI
            ROI=0D0
            PF=PHI(J)*P
            DO M1=1,MMAX2
               DL=DFLOAT(M1-1)*PF       
               CO=DCOS(DL)        
               A=2D0      
               IF(M1.EQ.1) A=1D0    
               ROI=ROI+A*RR(I,M1)*CO
            ENDDO
            IF (NSEP.EQ.1) THEN
               CALL MATR (MU(I),MU0,PF,P11)
               ROI=ROI+0.25D0*ALB*P11*MU0/(MU(I)+MU0)
            ENDIF
            RO(I,J)=ROI
         ENDDO
      ENDDO
      DO J=1,NPHI
         WRITE (14,1502) (RO(I,J), I=1,NMU)
      ENDDO
      STOP                                                              
      END                                                               

C********************************************************************   

C    CALCULATION OF THE PHASE FUNCTION FOR GIVEN MU, MU0, AND PHI-PHI0
C    AL1 - LEGENDRE EXPANSION COEFFICIENTS
C    MMAX1 - NUMBER OF COEFFICIENTS

      SUBROUTINE MATR (MU,MU0,PHI,P11)
      PARAMETER (LMAX1=700)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 AL1(LMAX1),MU,MU0
      COMMON /COEFF/ AL1,MMAX1
      PI=DACOS(-1D0)
      IF(DABS(PHI).LT.1D-7.OR.DABS(PHI-PI).LT.1D-7)
     &        PHI=PHI+1D-7
      IF (DABS(MU-MU0).LT.1D-10) MU0=MU0*0.999999999d0
      IF (DABS(MU+MU0).LT.1D-10) MU0=MU0*0.999999999d0
      SI=DSQRT(1D0-MU*MU)
      SI0=DSQRT(1D0-MU0*MU0)
      U=-MU*MU0+SI*SI0*DCOS(PHI)
      D6=DSQRT(6D0)*0.25D0
      F11=0D0
      P1=0D0
      PP1=1D0
      LMAX=MMAX1-1
      DO L1=1,MMAX1
         L=L1-1
         DL=DFLOAT(L)
         DL1=DFLOAT(L1)
         F11=F11+AL1(L1)*PP1
         IF(L.EQ.LMAX) GO TO 1
         PL1=DFLOAT(2*L+1)
         P=(PL1*U*PP1-DL*P1)/DL1
         P1=PP1
         PP1=P
      ENDDO   
    1 P11=F11
      RETURN
      END

C********************************************************************

      SUBROUTINE SPLINE (N,X,Y,B,C,D )                                  SPL00490
      IMPLICIT REAL*8 (A-H,O-Z)                                         SPL00500
      REAL*8 X(N),Y(N),B(N),C(N),D(N)                                   SPL00510
      NM1=N-1                                                           SPL00520
      IF (N.LT.2) RETURN                                                SPL00530
      IF (N.LT.3) GO TO 50                                              SPL00540
      D(1)=X(2)-X(1)                                                    SPL00570
      C(2)=(Y(2)-Y(1))/D(1)                                             SPL00580
      DO 10 I=2,NM1                                                     SPL00590
          D(I)=X(I+1)-X(I)                                              SPL00600
          DI=D(I)                                                       SPL00610
          B(I)=2D0*(D(I-1)+DI)                                          SPL00620
          C(I+1)=(Y(I+1)-Y(I))/DI                                       SPL00630
          C(I)=C(I+1)-C(I)                                              SPL00640
   10 CONTINUE                                                          SPL00650
      B(1)=-D(1)                                                        SPL00700
      NN1=N-1                                                           SPL00710
      B(N)=-D(NN1)                                                      SPL00720
      C(1)=0D0                                                          SPL00730
      C(N)=0D0                                                          SPL00740
      IF (N.EQ.3) GO TO 15                                              SPL00750
      C(1)= C(3)/(X(4)-X(2)) - C(2)/(X(3)-X(1))                         SPL00760
      NN2=N-2                                                           SPL00770
      NN3=N-3                                                           SPL00780
      C(N)= C(NN1)/(X(N)-X(NN2)) - C(NN2)/(X(NN1)-X(NN3))               SPL00790
      DD1=D(1)                                                          SPL00800
      C(1)=C(1)*DD1*DD1/(X(4)-X(1))                                     SPL00810
      DD1=D(NN1)                                                        SPL00820
      C(N) = -C(N)*DD1*DD1/(X(N)-X(NN3))                                SPL00830
   15 DO 20 I=2,N                                                       SPL00870
          II=I-1                                                        SPL00880
          T=D(II)/B(II)                                                 SPL00890
          B(I)=B(I)-T*D(II)                                             SPL00900
          C(I)=C(I)-T*C(II)                                             SPL00910
   20 CONTINUE                                                          SPL00920
      C(N)=C(N)/B(N)                                                    SPL00960
      DO 30 IB=1,NM1                                                    SPL00970
          I=N-IB                                                        SPL00980
          C(I)=(C(I)-D(I)*C(I+1))/B(I)                                  SPL00990
   30 CONTINUE                                                          SPL01000
      DDN=D(NM1)                                                        SPL01050
      B(N)= (Y(N)-Y(NM1))/DDN + DDN*(C(NM1)+2D0*C(N))                   SPL01060
      DO 40 I=1,NM1                                                     SPL01070
          II=I+1                                                        SPL01080
          CI=C(I)                                                       SPL01090
          CCI=C(II)                                                     SPL01100
          DI=D(I)                                                       SPL01110
          DDI=1D0/DI                                                    SPL01120
          B(I)= (Y(II)-Y(I))*DDI - DI*(CCI+2D0*CI)                      SPL01130
          D(I)=(CCI-CI)*DDI                                             SPL01140
          C(I)=3D0*CI                                                   SPL01150
   40 CONTINUE                                                          SPL01160
      C(N)=3D0*C(N)                                                     SPL01170
      D(N)=D(NN1)                                                       SPL01180
      RETURN                                                            SPL01190
   50 B(1)=(Y(2)-Y(1))/(X(2)-X(1))                                      SPL01200
      C(1)=0D0                                                          SPL01210
      D(1)=0D0                                                          SPL01220
      B(2)=B(1)                                                         SPL01230
      C(2)=0D0                                                          SPL01240
      D(2)=0D0                                                          SPL01250
      RETURN                                                            SPL01260
      END                                                               SPL01270

C********************************************************************

      DOUBLE PRECISION FUNCTION SEVAL ( N,U,X,Y,B,C,D )                 SEV00430
      IMPLICIT REAL*8 (A-H,O-Z)                                         SEV00450
      REAL*8 X(N),Y(N),B(N),C(N),D(N)                                   SEV00460
      DATA I/1/                                                         SEV00470
      IF (I.GE.N) I=1                                                   SEV00480
      IF (U.LT.X(I)) GO TO 10                                           SEV00490
      IF (U.LE.X(I+1)) GO TO 30                                         SEV00500
   10 I=1                                                               SEV00540
      J=N+1                                                             SEV00550
   20 K=(I+J)*0.5D0                                                     SEV00560
      XK=X(K)                                                           SEV00570
      IF (U.LT.XK) J=K                                                  SEV00580
      IF (U.GE.XK) I=K                                                  SEV00590
      IF (J.GT.I+1) GO TO 20                                            SEV00600
   30 DX=U-X(I)                                                         SEV00620
      SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))                            SEV00630
      RETURN                                                            SEV00640
      END                                                               SEV00650
Comments: