dtran_brdf.f

Go to the documentation of this file.

 
      subroutine 
     & diff_tran_corr(iphase,sun,view,
     &                  delphi,chlor,ta865,correct)
 
C     Iphase  is the aerosol model index (1=O50, 2=O70, ... 16=T99)
C     Iwave   is the wavelength index (1=412nm, 2=443nm, ... 6=670nm)
C     Sun     is the solar zenith angle (DEGREES)
C     View    is the satellite viewing angle (DEGREES)
C     Delphi  is the relative azimuth angle (RADIANS)  <---- NOTE NOTE NOTE
C     Chlor   is the chlorophyll concentration (mg/m^3)
C     Ta      is the aerosol optical depth for the GIVEN WAVELENGTH
C     Tr      is the Rayleigh optical depth for the GIVEN WAVELENGTH
C     Correct is vhe value of (t-t*)/t* for the GIVEN WAVELENGTH
 
 
c     The procedure is to apply the standard SeaWiFS/MODIS atmospheric correction 
c     algorithm, which provides tvLw.  During this procedure, the atmospheric 
c     correction algorithm yields two aerosol models (and the associated aerosol 
c     optical depths) and an interpolation parameter (Ratio?) that provides the 
c     fraction of each of the aerosol models to be used in the computation of 
c     the aerosol radiance.  Using the approximation tv = tv* and ts = ts* for 
c     each aerosol model (with the t*'s computed from the aerosol models and 
c     optical thicknesses and already provided in the Gordon and Wang algorithm) 
c     and interpolating between them using the same parameter (Ratio), the 
c     normalized water-leaving radiance and the chlorophyll concentration (Chl) 
c     is computed.  Then, using Chl, the subroutine above diff_tran_corr is 
c     entered for each visible wavelength (412 to 670 nm, SeaWiFS Bands) using 
c     the appropriate optical thicknesses for the each wavelength and the 
c     sun-sensor viewing geometry. This provides the diffuse transmittance 
c     correction factor Correct =(tv-tv*)/tv*.  The corrected diffuse 
c     transmittance is then 
c
c                       (Correct + 1) _ tv*.
c
c
c     This value of tv  is then used with the original tvLw to provide a better 
c     estimate of Lw.  It is expected that using the approximation tv = tv* will 
c     not have much influence on Chl so this probably does not need to be 
c     iterated.  If iteration is required, simply recomputed Chl and enter 
c     diff_tran_corr again.  This procedure will result in more precise 
c     values of Lw. 
 
 
 
      IMPLICIT NONE
      SAVE
 
      integer NRAD,NPHI,NWAVE,NPHASE,NGAUS,NNG,NG,NUM,NBIG
      integer I, J, M, Ibig, Iwave, Iphase
 
      parameter(nbig=10,nrad=31,nphi=4)       ! Note NRAD=31 NOT 91
      parameter(nphase=16,nwave=6)
 
      parameter(ngaus=2*nrad-1, nng=50, ng=2*nng)
      parameter(num=75)
 
      REAL Sun, View, Delphi, Ta865, Ta, Tr(NWAVE)
      REAL Correct(NWAVE), Chlor
      REAL  MBRDF(NWAVE,NRAD,NPHI)   ! Same as RAD1 in subroutine, 
                                     ! but has wave index RAD1(NWAVE,NRAD,NPHI)
      REAL  SmMBRDF(NWAVE,NBIG,NPHI) ! Shortened version with Mu averaging
      REAL PHSA(Nphase,Nwave,Nbig,NGAUS,NPHI), PHSR(Nbig,NGAUS,NPHI)
 
      REAL APHSRADA(NPHASE,NWAVE,NGAUS,NGAUS,NPHI)
 
      REAL MU(NGAUS),PDIV(NG),PWT(NG), THETA(NGAUS)
      REAL PHSRADA(NGAUS,NGAUS,NPHI), PHSRADR(NGAUS,NGAUS,NPHI)
      REAL TAUA_RAT(NPHASE,NWAVE)
 
      REAL  RAD1(NRAD,NPHI)
      REAL  tst(2),tdf(2)
 
      REAL TSTARR(NWAVE,NRAD), TSTARA(NPHASE,NWAVE,NRAD)
 
      REAL ya1,ya2,yr1,yr2
      REAL za1,za2,zr1,zr2
      REAL x1,z1,v1,w1
      REAL x2,z2,v2,w2
      REAL f1,f2,yl,order,fac,rfres
      REAL aint_p_a, aint_p_r, aint_pl_a, aint_pl_r
      REAL pterm_a, pterm_r, plterm_a, plterm_r
      REAL taur, taua, adelphi
      REAL t_star_r, t_star_a, t_diff_r, t_diff_a
      REAL tstar1, tstar2, tdiff1, tdiff2
      REAL ans1, ans2, ans, slope_ans, slope_tst, slope_tdf
      REAL PI, aindex, fresref, rad, ang, x 
 
      REAL tstartest1, tstartest2
 
      INTEGER JP, JDN, JUP, MAXPHI, MGAUS, IVIEW
 
      CHARACTER INFL1(2)*80,INFL2*80, DUMMY*80
 
c      common /comfour/ aphsrada, taua_rat
      common /comphase/ phsa, phsr, taua_rat
      COMMON /tstar/ tstarr, tstara 
 
      DATA tr /0.3132, 0.2336, 0.1547, 0.1330, 0.0957, 0.0446/
 
      rad(x)=x*pi/180.
      ang(x)=x*180./pi
 
      pi=4.0*atan(1.0)
 
      aindex = 1.334
 
      mgaus  = ngaus
      maxphi = nphi
 
      DO i = 1, ngaus   !NRAD
        theta(i)    = float(i-1)*90./float(nrad-1)
        mu(i)       = cos(rad(theta(i)))
      ENDDO
 
C     Load the appropriate aerosol quantities 
      call read_partial_phase_integrations ! After the first call to 
                                           ! this it just RETURNs
      call read_tstar              ! After the first call to this it just RETURNs
 
c       Load the Radiance/BRDF tables
 
      call morel_brdf(sun, chlor, mbrdf, smmbrdf)
 
C     Find where they are in the tables.
 
      DO i = 1, nrad
        IF( view .LT. theta(i) ) THEN
            iview = i
            GO TO 10
        ENDIF
      ENDDO
10    CONTINUE   ! USE IVIEW and IVIEW -1
 
 
399    FORMAT(a)
400    FORMAT( 5(3x, e12.6))
 
 
      ta = ta865
      DO iwave = 1, nwave
 
      DO jup = iview-1,iview
 
       adelphi = ang(delphi)
       jdn = mgaus + 1 - jup
 
C      First do the integrals that result from the scattering of Lw toward 
C      the sensor 
 
   
      ya1 = 0
c      YA2 = 0.
      yr1 = 0.
c      YR2 = 0.
      DO m = 1, maxphi
      order = float(m-1)
      x1 = 0.
c      Z1 = 0.
      v1 = 0.
c      W1 = 0.
       DO i = 1, nbig
        x2 = phsa(iphase,iwave,i,jup,m)*smmbrdf(iwave,i,m)
c        Z2 = PHSA(Iphase,Iwave,I,JUP,M)
        v2 = phsr(i,jup,m)*smmbrdf(iwave,i,m)
c        W2 = PHSR(I,JUP,M)
        x1 = x1 + x2
c        Z1 = Z1 + Z2
        v1 = v1 + v2
c       W1 = W1 + W2
       ENDDO
 
       IF(m.GT.1) THEN 
        x1 = x1*2
        v1 = v1*2
c       Z1 = 0.
c       W1 = 0.
       ENDIF
        ya1 = ya1 + cos(order*delphi)*x1  ! Index 1 has Rad1    (for t)
c        YA2 = YA2 + COS(ORDER*DELPHI)*Z1  ! Index 2 is w/o Rad1 (for t*)
        yr1 = yr1 + cos(order*delphi)*v1  
c        YR2 = YR2 + COS(ORDER*DELPHI)*W1
      ENDDO
 
C     Now do the integrals that result from the scattering of Lw back 
C     toward the surface with subsequint Fresnel reflection from the water. 
 
      za1 = 0
c      ZA2 = 0.
      zr1 = 0.
c      ZR2 = 0.
      DO m = 1, maxphi
      order = float(m-1)
      x1 = 0.
c      Z1 = 0.
      v1 = 0.
c      W1 = 0.
       DO i = 1, nbig
        x2 = phsa(iphase,iwave,i,jdn,m)*smmbrdf(iwave,i,m)
c        Z2 = PHSA(Iphase,Iwave,I,JDN,M)
        v2 = phsr(i,jdn,m)*smmbrdf(iwave,i,m)
c        W2 = PHSR(I,JDN,M)
        x1 = x1 + x2
c        Z1 = Z1 + Z2
        v1 = v1 + v2
c       W1 = W1 + W2
       ENDDO
       IF(m.GT.1) THEN 
        x1 = x1*2
        v1 = v1*2
c       Z1 = 0.
c       W1 = 0.
       ENDIF
        za1 = za1 + cos(order*delphi)*x1  ! Aerosol w    Rad1 for t
c        ZA2 = ZA2 + COS(ORDER*DELPHI)*Z1  ! Aerosol w/o  Rad1 for t*
        zr1 = zr1 + cos(order*delphi)*v1  ! Rayleigh w   Rad1 for t
c        ZR2 = ZR2 + COS(ORDER*DELPHI)*W1  ! Rayleigh w/o Rad1 for t*
      ENDDO
 
C     Compute the radiance in the viewing direction from Fourier coefficients
 
      yl = 0.
      DO m = 1, maxphi
        order = float(m-1)
        fac = 1.
        IF (m .GT. 1) fac = 2.
        yl = yl + mbrdf(iwave,jup,m)*fac*cos(order*delphi)
      ENDDO
 
C     Combine the two contributions above the integrals to form the diffuse 
C     reflectance in single scattering
 
      rfres = fresref(mu(jup), aindex)
 
c      aint_p_a  = (YA2 + rfres*ZA2)         /(1.-rfres)
      aint_pl_a = (ya1 + rfres*za1)/yl      /(1.-rfres)
c      aint_p_r  = (YR2 + rfres*ZR2)         /(1.-rfres)
      aint_pl_r = (yr1 + rfres*zr1)/yl      /(1.-rfres)
 
c      pterm_a  = (1.-aint_p_a/2. )/MU(JUP)     ! aint_p should be *omega0
      plterm_a = (1.-aint_pl_a/2.)/mu(jup)     !    "
c      pterm_r  = (1.-aint_p_r/2. )/MU(JUP)     
      plterm_r = (1.-aint_pl_r/2.)/mu(jup)      
 
c       t_star_a = exp(-pterm_a*ta)
        t_diff_a = exp(-plterm_a*ta*taua_rat(iphase,iwave))
c       t_star_r = exp(-pterm_r*tr)
        t_diff_r = exp(-plterm_r*tr(iwave))
        
            
        if (jup .EQ. iview-1) THEN 
c            tstar1 = t_star_a * t_star_r
            tstar1 = tstarr(iwave,jup)
     & * (tstara(iphase,iwave,jup)**(ta*taua_rat(iphase,iwave)))
            tdiff1 = t_diff_a * t_diff_r
            ans1   = (tdiff1 - tstar1)/(tstar1)
            
        else
c            tstar2 = t_star_a * t_star_r
            tstar2 = tstarr(iwave,jup)
     & * (tstara(iphase,iwave,jup)**(ta*taua_rat(iphase,iwave)))
            tdiff2 = t_diff_a * t_diff_r
            ans2   = (tdiff2 - tstar2)/(tstar2)
 
        endif
 
        ENDDO ! JUP
 
        slope_ans   = (ans2  -  ans1)/(theta(iview)-theta(iview-1))
        slope_tst   = (tstar2-tstar1)/(theta(iview)-theta(iview-1))
        slope_tdf   = (tdiff2-tdiff1)/(theta(iview)-theta(iview-1))
        correct(iwave)   =  ans1   + slope_ans*(view-theta(iview-1))
        
        ENDDO ! IWAVE
 
C        tst(iphase) = tstar1 + slope_tst*(view-theta(iview-1))
C        tdf(iphase) = tdiff1 + slope_tdf*(view-theta(iview-1))
 
C        ans = correct
 
c        print*, tstar1, tstar2
c        print*, adelphi, tst(Iphase), tdf(Iphase), ans 
 
      END
 
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
        function fresref(muair,  index)
 
        real muair, index
 
        theta1  = acos(muair)
        stheta1 = sin(theta1)
        stheta2 = stheta1/index
        theta2  = asin(stheta2)
 
        if(theta1 .gt. 0.) then 
 
                fresref = 0.5*(
     &         ( tan(theta1-theta2)/tan(theta1+theta2) )**2
     &        +( sin(theta1-theta2)/sin(theta1+theta2) )**2
     &             )
        else
                fresref = ( (index-1)/(index+1) )**2
        endif
        return
        end
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
 
C     *********************************************************
 
      subroutine morel_brdf(Sun, Chlor, MBRDF, SmMBRDF)
 
      implicit none
 
      integer NRAD,NPHI,NSUN,NCHL,NWAVE,NCASE,NBIG
      integer IRAD,IPHI,ISUN,ICHL,IWAVE,ICASE,IBIG
      integer I, M
 
      parameter(nbig=10,nrad=31,nphi=4)    ! Note NRAD=31 NOT 91
      parameter(nsun=6,nchl=6,nwave=6,ncase=nsun*nchl*nwave)
 
      REAL  BRDF(NWAVE,NSUN,NCHL,NRAD,NPHI)
      REAL  SmBRDF(NWAVE,NSUN,NCHL,NBIG,NPHI)
      REAL  Theta0(NSUN), Chl(NCHL), Wave(NWAVE), Thetav(NRAD)
      REAL  LChl(NCHL)
      REAL  SUN, Chlor, LChlor, chl_interp, sun_interp
      REAL  INTERP1, INTERP2
      REAL  MBRDF(NWAVE,NRAD,NPHI), SmMBRDF(NWAVE,NBIG,NPHI)
 
      CHARACTER*80  INFL_FOURIER31
      CHARACTER*80  INFL_FOURIER10
      CHARACTER*20  DUMMY
 
      CHARACTER filedir*255
      INTEGER   len, lenstr
 
      INTEGER IONCE
 
      SAVE
 
      DATA ionce /0/
 
      IF (ionce .EQ. 0) THEN
 
        call getenv('OCDATAROOT',filedir)
        if (filedir .eq. '') then
            write(*,*)
     .      '-E- : Environment variable OCDATAROOT undefined'
            call exit(1)
        endif
        len = lenstr(filedir)
 
      infl_fourier31 = filedir(1:len)//
     .  '/eval/common/dtran_brdf/NEW_Morel_NRAD31-1-EDITED'
      infl_fourier10 = filedir(1:len)//
     .  '/eval/common/dtran_brdf/NEW_Morel_SMALL-1-EDITED'
 
      OPEN(unit=11,file=infl_fourier31,status='UNKNOWN')
      OPEN(unit=12,file=infl_fourier10,status='UNKNOWN')
 
1     format(a)
2     format(' ', 13(2x,f10.6))
400   FORMAT( 5(3x, e12.6))
 
      DO iwave = 1, nwave
      DO isun  = 1, nsun
      DO ichl  = 1, nchl
      DO iphi = 1, nphi
         DO irad = 1, nrad
          brdf(iwave,isun,ichl,irad,iphi)=0.
         ENDDO
         DO ibig = 1, nbig
          smbrdf(iwave,isun,ichl,ibig,iphi)=0.
         ENDDO
      ENDDO
      ENDDO
      ENDDO
      ENDDO
 
C     READ IN THE FILES
 
      DO iwave = 1, nwave
      DO isun  = 1, nsun
      DO ichl  = 1, nchl
        Read(11, *) wave(iwave), theta0(isun), chl(ichl)
c        Write(6,*)  wave(iwave), theta0(isun), chl(ichl)
        DO iphi = 1, nphi
         Read(11, 1)   dummy
c         Write(6, 1)   DUMMY
         Read(11,400) (brdf(iwave,isun,ichl,irad,iphi),irad =1,nrad)
c         Write(6,400) (BRDF(iWAVE,iSUN,iCHL,iRAD,iPHI),irad =1,nrad)
        ENDDO
c        Read(12, *) wave(iwave), theta0(isun), chl(ichl)
        Read(12, *) wave(iwave), theta0(isun), chl(ichl)
c        Write(6,*)  wave(iwave), theta0(isun), chl(ichl)
        DO iphi = 1, nphi
         Read(12, 1)   dummy
c         Write(6, 1)   DUMMY
         Read(12,400) (smbrdf(iwave,isun,ichl,ibig,iphi),ibig =1,nbig)
c         Write(6,400) (SMBRDF(iWAVE,iSUN,iCHL,iBIG,iPHI),ibig =1,nbig)
        ENDDO
      ENDDO
      ENDDO
      ENDDO
 
      ionce = 1
      ENDIF
 
      DO i = 1, nchl
       lchl(i) = alog10(chl(i))
      ENDDO
 
C     For sun angle interpolation
      DO i = 1, nsun
        IF( sun .LT. theta0(i) ) THEN
            isun = i
            GO TO 10
        ENDIF
      ENDDO
10    CONTINUE   ! USE ISUN and ISUN -1
 
C     For Chl interpolation
      DO i = 1, nchl
        IF( chlor .LT. chl(i) ) THEN
            ichl = i
            GO TO 11
        ENDIF
        ichl = nchl
      ENDDO
11    CONTINUE   ! USE ICHL and ICHL -1
 
c       Compute the log10 of the chlorophyll
        lchlor = alog10(chlor)
 
 
C       Interpolate in the tables to get the water radiance
 
        sun_interp = (sun   - theta0(isun-1) )/(theta0(isun)-theta0(isun-1))
        chl_interp = (lchlor - lchl(ichl-1) )/(lchl(ichl)-lchl(ichl-1))
c        PRINT*, 'Interpolation terms ', sun_interp, chl_interp
        
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
 
        DO iwave = 1, nwave
        DO i     = 1, nrad
        DO m     = 1, nphi
            interp1 = (1.-sun_interp)*(1.-chl_interp)*
     &                brdf(iwave,isun-1,ichl-1,i,m)
            interp1 = interp1 + sun_interp*(1.-chl_interp)*
     &                brdf(iwave,isun,ichl-1,i,m)
            interp1 = interp1 + (1.-sun_interp)*chl_interp*
     &                brdf(iwave,isun-1,ichl,i,m)
            interp1 = interp1 + sun_interp*chl_interp*
     &                brdf(iwave,isun,ichl,i,m)
            mbrdf(iwave, i, m) = interp1
 
            IF(i .LE. nbig) THEN
              interp2 = (1.-sun_interp)*(1.-chl_interp)*
     &                smbrdf(iwave,isun-1,ichl-1,i,m)
              interp2 = interp2 + sun_interp*(1.-chl_interp)*
     &                smbrdf(iwave,isun,ichl-1,i,m)
              interp2 = interp2 + (1.-sun_interp)*chl_interp*
     &                smbrdf(iwave,isun-1,ichl,i,m)
              interp2 = interp2 + sun_interp*chl_interp*
     &                smbrdf(iwave,isun,ichl,i,m)
              smmbrdf(iwave, i, m) = interp2
            ENDIF
c        print*,  'Iwave,I,M,MBRDF = ',Iwave,I,M,MBRDF(Iwave, I, M)
        ENDDO
        ENDDO
        ENDDO
 
        DO iwave = 1, nwave
        DO i     = 1, nbig
        DO m     = 1, nphi
c        print*,  'Iwave,I,M,SmMBRDF = ',Iwave,I,M,SmMBRDF(Iwave, I, M)        
        ENDDO
        ENDDO
        ENDDO
 
 
      RETURN
      END 
      
 
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
      subroutine read_partial_phase_integrations
 
      IMPLICIT NONE
 
      integer NRAD,NPHI,NWAVE,NPHASE,NGAUS,NBIG
 
      parameter(nrad=31,nphi=4,nbig=10,ngaus=2*nrad-1)  !Note NRAD=31 NOT 91
      parameter(nphase=16,nwave=6)
 
      integer I, J, JUP, M, Iwave, Iphase, Ibig, IONCE
 
      REAL PHSA(Nphase,Nwave,Nbig,NGAUS,NPHI), PHSR(Nbig,NGAUS,NPHI)
      REAL TAUA_RAT(NPHASE,NWAVE)
 
      CHARACTER INFL_AER*512,INFL_RAY*512, DUMMY*80
      CHARACTER INFL_RAT*512
 
      character filedir*255
      integer   len, lenstr
 
      common /comphase/ phsa, phsr, taua_rat
 
      DATA  ionce /0/
 
      Save
 
      IF (ionce .EQ. 0) THEN
 
         call getenv('OCDATAROOT',filedir)
         if (filedir .eq. '') then
            write(*,*)
     .      '-E- : Environment variable OCDATAROOT undefined'
            call exit(1)
         endif
         len = lenstr(filedir)
         filedir = filedir(1:len)//'/eval/common/dtran_brdf/'
         len = lenstr(filedir)
 
      infl_aer = filedir(1:len)//'Aerosols_Partial_Inegr.dat'
      infl_ray = filedir(1:len)//'Rayleigh_Partial_Inegr.dat'
      infl_rat = filedir(1:len)//'spec_var_EDITED.dat'
 
      OPEN(unit=11,file=infl_aer,status='OLD')
      OPEN(unit=12,file=infl_ray,status='OLD')
      OPEN(unit=15,file=infl_rat,status='OLD')
 
399    FORMAT(a)
400    FORMAT( 5(3x, e12.6))
 
 
      DO iphase = 1,nphase
       DO iwave = 1,nwave
        DO m = 1, nphi
         DO jup = 1, ngaus  ! Read in for thetav = 0, 3, 6, 9, etc.
           Read(11,399) dummy
c           Write(6,399) Dummy
           Read(11,400) (phsa(iphase,iwave,ibig,jup,m), ibig = 1, nbig)
c           Write(6,400) (PHSA(Iphase,Iwave,Ibig,JUP,M), Ibig = 1, NBIG)
         ENDDO      !JUP
        ENDDO       !M
       ENDDO        !Iwave
      ENDDO         !Iphase
 
      CLOSE(11)
 
       DO m = 1, nphi
        DO jup = 1, ngaus  ! Read in for thetav = 0, 3, 6, 9, etc.
           Read (12,399) dummy
c           Write( 6,399) Dummy
           Read (12,400) (phsr(ibig,jup,m), ibig = 1, nbig)
c           Write( 6,400) (PHSR(Ibig,JUP,M), Ibig = 1, NBIG)
        ENDDO      !JUP
       ENDDO       !M
 
       CLOSE(12)
 
C     READ THE TAUA RATIOS (Lambda to 865)
 
        do iphase = 1, nphase
         read(15,399) dummy
C         write(6,399) Dummy
         do iwave = 1, nwave
            read(15,*) j, taua_rat(iphase,iwave)
C            write(6,*) iphase, iwave, j, taua_rat(iphase,iwave)
         enddo
        enddo
        close(15)
 
       ionce = 1
       ENDIF
      RETURN
      END
 
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
      SUBROUTINE read_tstar
 
      IMPLICIT NONE
 
      integer NRAD,NPHI,NWAVE,NPHASE
      integer I, J, Iwave, Iphase, IONCE
 
      parameter(nrad=31)       ! Note NRAD=31 NOT 91
      parameter(nphase=16,nwave=6)
 
      REAL TSTARR(NWAVE,NRAD), TSTARA(NPHASE,NWAVE,NRAD)
 
      CHARACTER*80 Dummy
 
      character filedir*255
      integer   len, lenstr
 
      COMMON /tstar/ tstarr, tstara 
 
      DATA ionce /0/
 
      SAVE
 
 397   Format(' ' , '  Iwave = ', i2)
 398   Format(' ' , ' Iphase = ', i2, '  Iwave = ', i2)
 399   FORMAT(a)
 400   FORMAT( 5(3x, e12.6))
      
 
      If (ionce .EQ. 0) THEN
 
         call getenv('OCDATAROOT',filedir)
         if (filedir .eq. '') then
            write(*,*)
     .      '-E- : Environment variable OCDATAROOT undefined'
            call exit(1)
         endif
         len = lenstr(filedir)
         filedir = filedir(1:len)//'/eval/common/dtran_brdf/'
         len = lenstr(filedir)
 
       OPEN(unit=21,file=filedir(1:len)//'tstar_rayleigh.dat',
     .      status='UNKNOWN')
       OPEN(unit=22,file=filedir(1:len)//'tstar_aerosol.dat',
     .      status='UNKNOWN')
 
        DO iphase = 1, nphase
         DO iwave = 1, nwave
           If(iphase .eq.1) then 
              Read(21,399) dummy
              Read(21,400) (tstarr(iwave, j), j = 1, nrad-3)
c              Write(6,397) Iwave
c              Write(6,400) (TSTARR(Iwave, J), J = 1, Nrad-3)
           Endif
              Read(22,399) dummy
              Read(22,400) (tstara(iphase,iwave, j), j = 1, nrad-3) 
c              Write(6,398) Iphase, Iwave
c              Write(6,400) (TSTARA(Iphase,Iwave, J), J = 1, Nrad-3)
         ENDDO
        ENDDO   
 
        Close(21)
        Close(22)
        
        ionce = 1
 
        RETURN
      ELSE
        RETURN
      ENDIF
 
      END
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC