!------------------------------------------------------------
!
!  MODULE OSS_ADDBL: contains items needed for the 
!                OSS Scattering Adding-Doubling Module
!
!--------------------------------------------------------------
 MODULE oss_addbl
  PRIVATE 
  
  public:: scattRT
  
 !------------------------------------------------------------------------
 !    Profile Parameters (should match with oss_ir_module)
 !------------------------------------------------------------------------
      integer    :: mxlay,maxcld
      parameter (mxlay=100,maxcld=10)
 !------------------------------------------------------------------------
 !    Multiple Scattering Parameters
 !------------------------------------------------------------------------
        real accur
        integer mxcang,maxcmu,mxang,mxang2
        parameter(mxcang=64)
        parameter (maxcmu=2*mxcang)
        parameter(mxang=mxcang+1,mxang2=mxang*mxang)

        data accur/1.e-04/
        
CONTAINS
  
  SUBROUTINE scattRT(tautot,xG,btemp,tabsHydr,tscatt,gasym,emis,fBeam,vn,nLay,nObs,&
       sBeam,nStr,umu,umu0,lookup,delphi,iflux,flxNetMono, &
       rad,xkt,xkEmRf)
        
    !---Wrapper to interface oss_module with ossScat.  Some calling arguments
    !      are placeholders.
    USE constants, only : MISSING_REAL
    
    !--Input/output variables
    REAL                             :: vn,emis,fBeam,umu,umu0,delphi,rad
    INTEGER                          :: nLay,nObs,nStr,iflux
    LOGICAL                          :: lookup,sbeam
    INTEGER,DIMENSION(:),ALLOCATABLE :: lTyp
    REAL, DIMENSION(:)               :: xG,tautot,flxNetMono
    REAL, DIMENSION(:)               :: tabsHydr,tscatt,gasym
    REAL, DIMENSION(:)               :: xkt
    REAL                             :: xkEmRf
    !---Internal variables
    REAL*8                 :: rvn
    REAL                   :: tau(mxlay),cldabs
    REAL                   :: albedo,btemp,gg,hl(0:maxcmu),pmom(0:maxcmu,mxlay)
    REAL, PARAMETER        :: scatLB=1e-8
    INTEGER                :: l,k,ik,ikm1
    LOGICAL                :: lamber,plank,deltam
    REAL,DIMENSION(:),ALLOCATABLE :: tscattIn
    
    !---Variable intialization
    data lamber/.true./,deltam/.true./
    plank=.false.
    if(vn.lt.5000.)plank=.true.
        
    if(.not.sbeam)fBeam=0.
    
    albedo=1.-emis
    rvn = vn
    
    xkt  = MISSING_REAL
    xkem = MISSING_REAL
    
    allocate(lTyp(nLay),tscattIn(size(tscatt)))
    tscattIn = 0.
    pmom     = 0.
    lTyp = 0
    !---Compute optical properties of default layers:
    do l=1,nLay
       !---Absorption optical depth (water vapor, fixed gases and hydrometeor):
       tau(l) =tautot(l) + tabsHydr(l)
       if (tscatt(l).gt.scatLB) lTyp(l)=1
       if (lTyp(l).gt.0.) then
          !---Phase function:
          pmom(0,l)=1.
          gg=gasym(l)
          do k=1,nstr
             pmom(k,l)=pmom(k-1,l)*gg
          end do
       end if
    end do
    
    tscattIn = tscatt
    call ossscat( nlay, tau, tscattIn, pmom, xG, lTyp, & 
         rvn, nstr, nobs, umu, fBeam, umu0, delphi, lamber, albedo, hl, &
         btemp, deltam, plank, sbeam, lookup, rad)
    
    deallocate(lTyp,tscattIn)
    
    RETURN  
    
  END SUBROUTINE SCATTRT
      
      SUBROUTINE OSSSCAT( NLAY, TAUABS, TAUSCAT, ALPHA, TEMP, LTYPIN, &
         VN, NSTR, LOBS, UMU, FBEAM, UMU0, DELPHI, LAMBER, ALBEDO, HL, &
         TSFC, DELTAM, PLANKIN, SBEAMIN, LOOKUP, RADOUT)

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! **** ***********   VERSION 2.0 :    21 MARCH 1995   *****************
!
!                       ALGORITHM:    J.L. MONCET
!                                     S.A. CLOUGH
!
!                  IMPLEMENTATION:    J.L. MONCET
!                                     D.A. PORTMAN
!
!               ATMOSPHERIC AND ENVIRONMENTAL RESEARCH INC.
!               840 MEMORIAL DRIVE,  CAMBRIDGE, MA   02139
!
!----------------------------------------------------------------------
!
!               WORK SUPPORTED BY:    THE ARM PROGRAM
!                                     OFFICE OF ENERGY RESEARCH
!                                     DEPARTMENT OF ENERGY
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL*8 VN
      REAL TAUABS(MXLAY),TAUSCAT(MXLAY)
      REAL ALPHA(0:2*MXCANG,MXLAY)
      REAL TEMP(0:MXLAY)
!---- REAL UMU,UMU0,DELPHI,FBEAM,ACCUR
      REAL UMU,UMU0,DELPHI,FBEAM
      REAL ALBEDO,HL(1),TSFC,RADOUT
      INTEGER NLAY,NSTR,LTYPIN(MXLAY),LTYP(MXLAY),LOBS
      LOGICAL DELTAM,PLANKIN,PLANK,SBEAMIN,SBEAM,LOOKUP,LAMBER
!
      REAL TAUEXT(MXLAY),SSALB(MXLAY),F0(0:MXLAY)
      real alpha1(0:2*mxcang),rad(mxang)
      LOGICAL REFLB,REFLT
!
      REAL TTOP(1),TBOT(1)

      COMMON /LAYDAT/ RLAY(MXANG2),TLAY(MXANG2), EMPL(MXANG),EMML(MXANG)
      COMMON /ACMDWN/ RTOP(MXANG2),EMTOP(MXANG)
      COMMON /ACMUP / RBOT(MXANG2),EMBOT(MXANG)

      COMMON /ANGL  / PTH(MXANG),GMU(MXANG),GWT(MXANG),GMU0,&
                     YLM(0:2*MXCANG-1,MXANG),YLM0(0:2*MXCANG-1),&
                     GMU0VEC(1)

      COMMON /LPAR  / PLANK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /CONSTN/ PI,CL2INV

!
      DATA INIT/1/
      SAVE INIT
!
      PLANK=PLANKIN
      SBEAM=SBEAMIN
      LTYP(1:NLAY)=LTYPIN(1:NLAY)
!
      IF(INIT.EQ.1)THEN
! ======================================================================
! === INITIALIZATION ===================================================
! ======================================================================
         PI=ACOS(-1.)
         CL2INV=1./LOG(2.)
!
         NCANG=NSTR/2
!
!-----Generate discrete angles and associated weights for
!-----Double-Gauss" quadrature (see Stamnes, 1981)
!
         CALL QGAUSN(NCANG,GMU,GWT)
!
         NANG=NCANG+1
         N2ANG=2*NANG
         NANG2=NANG*NANG
!
         DO 120 I=1,NCANG
  120    PTH(I)=1./GMU(I)
!
         INIT=0
      END IF
!
! *** Insert user angle.
!
      GMU(NANG)=UMU
      PTH(NANG)=1./UMU
!
! *** Set NAZ
!
      NAZ=0
      IF(SBEAM)THEN
         GMU0=UMU0
         IF((1.-GMU0).GT.1.E-05.AND. (1.-GMU(NANG)).GT.1.E-05)NAZ=N2ANG-3
      END IF
!
! *** Set control parameters for CHARTS run
!
      IF(LAMBER) THEN
         ISRFL=-1
         IF(ALBEDO.LE.1.E-6)ISRFL=0
      ELSE
         STOP 'OSSSCAT: This version only supports Lambertian surfaces'
      ENDIF
!
      CALL SETCTL(NLAY,LTYP,LOBS,SBEAM,PLANK,ISRFL,LOOKUP,NAZ,&
         NBOT0,NBOT,NTOP0,NTOP)
!
      N1=NTOP0
      N2=NBOT0
!
! *** Delta-M scaling
!
      F=0.
      DO I=1,NLAY
         IF(DELTAM)THEN
            F=ALPHA(N2ANG-2,I)
            TAUSCAT(I)=TAUSCAT(I)*(1.-F)
         END IF
!
         DO 50 L=0,N2ANG-3
            ALPHA(L,I)=(2*L+1)*(ALPHA(L,I)-F)/(1.-F)
 50      CONTINUE
      ENDDO
!
      DO 330 L=1,NLAY
         TAUEXT(L)=TAUABS(L)
         IF(LTYP(L).GT.0)THEN
            TAUEXT(L)=TAUEXT(L)+TAUSCAT(L)
            IF(TAUEXT(L).LE.0.)THEN
               SSALB(L)=0.
            ELSE
               SSALB(L)=TAUSCAT(L)/TAUEXT(L)
            END IF
         END IF
  330 CONTINUE
!
      IF(SBEAM)THEN
!
!-----Compute solar beam extinction (current version does not treat
!-----specularly reflected beam)
!
         SUMEX=0.
         F0(0)=FBEAM
         DO 350 L=1,NLAY
!-----Add scattering component to total extinction
            SUMEX=SUMEX+TAUEXT(L)
!-----Calculate solar irradiance at top of scattering layer
            F0(L)=FBEAM*EXP(-SUMEX/GMU0)
  350    CONTINUE
      END IF
!
! =====================================================================
! === BEGIN LOOP OVER AZIMUTHAL ANGLES ================================
! =====================================================================
      KCONV=0
      DO 700 MAZ=0,NAZ
       EMBOT=0.
       EMTOP=0.
       RTOP=0.
       RBOT=0.

!
! *** Generate normalized Legendre polynomials for computational angles
! *** (GMU) and for incident solar beam angle cosine (GMU0).
!
      CALL LEPOLY_OSS(NANG,MAZ,2*MXCANG-1,N2ANG-3,GMU,YLM)
!     --Put GMU0 in vector form to prevent compiler error:
      GMU0VEC(1)=GMU0
      IF(SBEAM)CALL LEPOLY_OSS(1,MAZ,2*MXCANG-1,N2ANG-3,GMU0VEC,YLM0)
!      IF(SBEAM)CALL LEPOLY_OSS(1,MAZ,2*MXCANG-1,N2ANG-3,GMU0,YLM0)
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     LOOP OVER LAYERS THE ABOVE OBSERVER LEVEL - MERGE DOWN CASE      +
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      REFLT=.FALSE.
!
      DO 500 LNUM=N1,LOBS-1
!
         IF(LTYP(LNUM).LE.0)THEN
!-----------------------------------------------------------------------
!     CLEAR LAYER CASE: 
!-----------------------------------------------------------------------
! *** Accumulate partial results for successive clear layers "LAYDAT".
! *** Final merge with results contained in "ACMDWN" done at the end,
! *** when (LTYP=-99 OR -100).
            CALL CLRLAY(TAUEXT(LNUM),TEMP(LNUM-1),&
               TEMP(LNUM),REFLT,LTYP(LNUM),VN)
            IF(LTYP(LNUM).LE.-99) CALL ADDCLR(RTOP,TTOP,EMTOP,REFLT,PLANK)
         ELSE
!-----------------------------------------------------------------------
!     CLOUDY LAYER CASE: 
!-----------------------------------------------------------------------
! *** Calculate reflection/transmission matrices (in RLAY and TLAY) and 
! *** internal sources (EMPL-upward- and EMML-downward) for current layer
!
            do k=0,n2ang-2
               alpha1(k)=alpha(k,lnum)
            end do
            CALL MSLAY(TAUEXT(LNUM),SSALB(LNUM),TEMP(LNUM-1),TEMP(LNUM),&
               F0(LNUM-1),RLAY,TLAY,EMPL,EMML,MAZ,VN,ALPHA1,0)
!
! *** Accumulate results in ACMWDN
!
            CALL ADDMS(RTOP,TTOP,EMTOP,REFLT)
         END IF
  500 CONTINUE
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     LOOP OVER LAYERS BELOW THE OBSERVER LEVEL - MERGE UP CASE        +
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! *** Initialize RBOT and EMBOT
!
      REFLB=.FALSE.
      CALL SURFCE(MAZ,ISRFL,TSFC,ALBEDO,F0(NLAY),RBOT,EMBOT,REFLB,VN)
!     
      DO 600 LNUM=N2,LOBS,-1
!     
         IF(LTYP(LNUM).LE.0)THEN
!-----------------------------------------------------------------------
!     CLEAR LAYER CASE: 
!-----------------------------------------------------------------------
! *** Accumulate partial results for successive clear layers "LAYDAT".
! *** Final merge with results contained in "ACMUP" done at the end,
! *** when (LTYP=-99 OR -100).
            CALL CLRLAY(TAUEXT(LNUM),TEMP(LNUM),TEMP(LNUM-1),REFLB,LTYP(LNUM),VN)
            IF(LTYP(LNUM).LE.-99) CALL ADDCLR(RBOT,TBOT,EMBOT,REFLB,PLANK)
         ELSE
!-----------------------------------------------------------------------
!     CLOUDY LAYER CASE: 
!-----------------------------------------------------------------------
! *** Calculate reflection/transmission matrices (in RLAY and TLAY) and 
! *** internal sources (EMPL-upward- and EMML-downward) for current layer
!
            do k=0,n2ang-2
               alpha1(k)=alpha(k,lnum)
            end do
            CALL MSLAY(TAUEXT(LNUM),SSALB(LNUM),TEMP(LNUM),TEMP(LNUM-1),&
               F0(LNUM-1),RLAY,TLAY,EMPL,EMML,MAZ,VN,ALPHA1,1)
!     
! *** Accumulate results.
!     
            CALL ADDMS(RBOT,TBOT,EMBOT,REFLB)
         END IF
  600 CONTINUE
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     DO FINAL MERGE  - Output radiances in RAD                        +
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      IF (LOOKUP) THEN
         CALL MRGRAD(RTOP,EMTOP,RBOT,EMBOT,REFLT,REFLB,RAD)
      ELSE
         CALL MRGRAD(RBOT,EMBOT,RTOP,EMTOP,REFLB,REFLT,RAD)
      END IF
! =====================================================================
! === END OF PANEL LOOP ===============================================
! =====================================================================
      IF(MAZ.EQ.0)THEN
         RADOUT=RAD(NANG)
! *** (Reset control parameters for MAZ >0)
         PLANK=.FALSE.
         N1=NBOT
         N2=NTOP
      ELSE
! *** Compute radiance scaling for azimuthal order MAZ
         AZTERM=RAD(NANG)*COS(MAZ*DELPHI*PI/180.)
         RADOUT=RADOUT+AZTERM
!
         IF(RADOUT.NE.0.)THEN
            AZERR=ABS(AZTERM)/ABS(RADOUT)
            IF(AZERR.LE.ACCUR)KCONV=KCONV+1
         END IF
!
         IF (KCONV.GE.2) GOTO 800
!            
      END IF
  700 CONTINUE
! =====================================================================
! === END OF LOOP OVER AZIMUTHAL ANGLES ===============================
! =====================================================================
  800 RETURN
!
      END subroutine ossscat
!
      SUBROUTINE SETCTL(NLAY,LTYP,LOBS,SBEAM,PLANK,ISRFL,&
         LOOKUP,NAZ,NBOT0,NBOT,NTOP0,NTOP)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     This subroutine adjusts LTYP for the clear layers and adjusts the
!     atmospheric boundaries to avoid in many cases treating systematically 
!     all the input layers in the radiative transfer calculations (for 
!     instance, in the solar case one can ignore the uppermost layers if 
!     they are non-reflecting since they do not produce a diffuse radiation 
!     field. Similarly, in presence of a Lambertian surface, there is no 
!     contribution from the surface and the clear layers adjacent to the 
!     surface for high order azimuthal components. CHARTS calculations 
!     will in this case, be restricted to the portion of the atmosphere 
!     between the top and bottom scattering layers). This subroutine also 
!     resets ISBEAM, IPLANK and NAZ in order to save on the amount of 
!     computation in some special cases (e.g. when the observer faces a 
!     non-reflecting atmosphere).
!
!     Values of LTYP for clear layers are modified to indicate the 
!     relative position of the each layer within a set of contiguous
!     clear layers above and below the observer level. LTYP values
!     are assigned values as follows:
!
!     LTYP = -1   - first clear layer
!          =  0   - middle layers 
!          = -99  - last clear layer
!          = -100 - single clear layer
!
!     The input values of LTYP for cloudy layers indicate the actual cloud
!     model used for each layer. LTYP is modified in this routine to provide 
!     a unique index value (>1) for each cloud layer.
!
!     Output variables:
!
!     NBOT0= lowest atmospheric layer below observer level at  which
!            one should start the radiative transfer calculations for
!            MAZ=0
!     NBOT = same as NBOT0 for MAZ >0
!     NTOP0= highest atmospheric layer above observer level at  which
!            one should start the radiative transfer calculations for
!            MAZ=0
!     NTOP = same as NTOP0 for MAZ >0
!
!     REVISED:  21 MARCH 1995
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      INTEGER LTYP(NLAY)
      LOGICAL PLANK,SBEAM,LOOKUP
      LOGICAL REFLS0,REFLS,REFLT,REFLB0,REFLB
      DATA REFLS0,REFLS/.FALSE.,.FALSE./
!
! *** Surface reflectance
!
      IF(ISRFL.NE.0)REFLS0=.TRUE.
      IF(ISRFL.GT.0)REFLS=.TRUE.
!
! *** Set LTYP for clear layers above observer level (MERGE DOWN case)
!     
      NTOP=1
      REFLT=.FALSE.
      DO 10 I=1,LOBS-1
         IF(LTYP(I).EQ.0)THEN
!
            IF(I.EQ.1)THEN
               LTYP(I)=-1
            ELSE IF(LTYP(I-1).GT.0)THEN
               LTYP(I)=-1
            END IF
!
            IF(I.EQ.LOBS-1)THEN
               LTYP(I)=-99+LTYP(I)
            ELSE IF(LTYP(I+1).GT.0)THEN
               LTYP(I)=-99+LTYP(I)
            END IF
!
            IF (.NOT.REFLT) NTOP=NTOP+1
         ELSE
            REFLT=.TRUE.
         END IF
  10  CONTINUE
!
! *** Set LTYP for clear layers below observer level (MERGE UP case)
!
      NBOT=NLAY
      REFLB=.FALSE.
      DO 20 I=NLAY,LOBS,-1
         IF(LTYP(I).EQ.0)THEN
!
            IF(I.EQ.NLAY)THEN
               LTYP(I)=-1
            ELSE IF(LTYP(I+1).GT.0)THEN
               LTYP(I)=-1
            END IF
!
            IF(I.EQ.LOBS)THEN
               LTYP(I)=-99+LTYP(I)
            ELSE IF(LTYP(I-1).GT.0)THEN
               LTYP(I)=-99+LTYP(I)
            END IF
!
            IF(.NOT.REFLB)NBOT=NBOT-1
         ELSE
            REFLB=.TRUE.
         END IF
  20  CONTINUE
!
      REFLB0=REFLB.OR.REFLS0
      REFLB=REFLB.OR.REFLS
!
! *** Adjust boundaries for general solar case
!
      NBOT0=NBOT
      NTOP0=NTOP
      IF(REFLS0)NBOT0=NLAY
      IF(REFLS)NBOT=NLAY
!
! *** Adjust boundaries for thermal case
!
      IF(PLANK)THEN
         NBOT0=NLAY
         NTOP0=1
      END IF
!
! *** Special case: No forward reflection - shrink atmospheric
! *** boundaries and force SBEAM = .FALSE. (assume that observer 
! *** never looks directly into the sun.)
!
      IF(LOOKUP.AND..NOT.REFLT)THEN
         SBEAM=.FALSE.
         NAZ=0
!         ISRFL=-99
         NBOT0=LOBS-1
      END IF
!
      IF(.NOT.LOOKUP)THEN
         IF(.NOT.REFLB)NAZ=0
         IF(.NOT.REFLB0)THEN
            SBEAM=.FALSE.
            NTOP0=LOBS
         END IF
      END IF
!
! *** Eliminate trivial case 
!
!      IF ((LOOKUP.AND.LOBS.EQ.1).OR.&
!         (.NOT.LOOKUP.AND..NOT.REFLB0.AND.LOBS.EQ.NLAY+1))&
!         STOP 'CASE NOT TREATED'
      RETURN
      END subroutine setctl
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     OSS (Optimal Spectral Sampling)
!     $Name:  $ 
!     $Id: oss_addbl.f90,v 1.4 2006/12/15 21:04:57 raschbre Exp $ 
!     Copyright AER, Inc., 2002, 2003. All rights Reserved.
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE SURFCE(MAZ,ISRFL,TSFC,SALB,F0,R,EM,REFLS,VN)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine initializes the reflection matrix (R) and emission 
!     vector (EM) for the lower atmosphere by applying boundary conditions
!     at the surface. The current version treats Lambertian surfaces only.
!
!     MODIFIED:  21 MARCH 1995 to implement spectral interpolation
!                of surface properties within panel
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL R(NANG,NANG),EM(NANG)
      REAL*8 VN
      LOGICAL PLNCK,SBEAM,REFLS

      COMMON /ANGL  / PTH(MXANG),GMU(MXANG),GWT(MXANG),GMU0,&
                     YLM(0:2*MXCANG-1,MXANG),YLM0(0:2*MXCANG-1)
      COMMON /CONSTN/ PI,CL2INV
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2

!
!
      REFLS=.FALSE.
!
      CALL VZERO(EM,NANG)
!
      IF(ISRFL.EQ.-99)RETURN
!      IF(ISRFL.GT.0)GOTO 100
!
! *** LAMBERTIAN SURFACE
!
      IF (MAZ.GT.0)RETURN
!     
      RS=0.
!
      IF(ISRFL.EQ.-1)THEN
!     
! *** Calculate surface reflectance across panel
!
         RS=SALB
         REFLS=.TRUE.
      ELSE IF (ISRFL.EQ.2)THEN
         STOP 'OSSCAT: Subroutine SURFCE - option not supported'
      END IF
!     
! *** Calculate reflection matrix
!
      IF(REFLS)THEN
         DO 40 J=1,NANG-1
            FAC=2.*GMU(J)*GWT(J)
            DO 40 I=1,NANG
  40           R(J,I)=RS*FAC

            DO 45 I=1,NANG
  45           R(NANG,I)=0.
      END IF
!
! *** Initialize EM by combining thermal emission and reflection of direct 
! *** solar beam at the surface
!
      IF(PLNCK)THEN
         B=WNPLAN(VN,TSFC)
         WK=(1.-RS)*B
!     
         DO 60 I=1,NANG
  60        EM(I)=WK
      END IF
!
      IF(SBEAM.AND.REFLS)THEN
         FAC=GMU0/PI
         WK=RS*F0*FAC
!
         DO 75 I=1,NANG
  75        EM(I)=EM(I)+WK
      END IF
!
! *** General case (not implemented)
!
  100 RETURN
      ENDsubroutine surfce
!
      SUBROUTINE CLRLAY(TAU,TA,TB,REFLG,INIT,VN)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the transmission and emission calculations
!     for a clear slab. The routine processes NLIM spectral points.
!     Transmission and forward and backward emission for each successive 
!     clear layer are accumulated respectively in TLAY, EMPL and EMML.
!     If no thermal emission the routine accumulates the optical depths 
!     in TLAY and calculates the transmission only when the last clear layer 
!     of a series is reached. Linear-in-tau approximation is used for the 
!     thermal source and it is assumed that the Planck function varies 
!     linearly across the spectral interval spanned by the NLIM spectral 
!     points.
!
!     REVISED:  11 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL*8 VN
      LOGICAL REFLG,PLNCK,SBEAM
!
      COMMON /ANGL  / PTH(MXANG),GMU(MXANG),GWT(MXANG),GMU0,&
                     YLM(0:2*MXCANG-1,MXANG),YLM0(0:2*MXCANG-1)
      COMMON /CONSTN/ PI,CL2INV
      COMMON /LAYDAT/ RLAY(MXANG2),TLAY(MXANG2),&
                     EMPL(MXANG),EMML(MXANG)
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2

      COMMON /LOCWK / TX(MXANG),EMP(MXANG),EMM(MXANG),WK(MXANG)
!
!     Scale optical depths for all incident angles including observer
!     angle.
!
      DO 20 I=1,NANG
         WK(I)=TAU*PTH(I)
  20  CONTINUE
!
! *** THERMAL CASE. 
! 
      IF(PLNCK)THEN
!     
!     Calculate vertically averaged Planck function for the layer and
!     correction term for linear-in-Tau approximation at each spectral 
!     point
!
         BA=WNPLAN(VN,TA)
         BB=WNPLAN(VN,TB)
         BAV=.5*(BB+BA)
         DELB=.5*(BB-BA)
! Or omit linear in tau
!        bav = wnplan(vn,(ta+tb)/2.)
!        delb = 0.
!
!     Compute transmittance and thermal emission
!
         DO 50 N=1,NANG
            TX(N)=EXP(-WK(N))
            EM0=1.-TX(N)
            IF(WK(N).GT.1.E-03)THEN
               CCI=2./WK(N)
               EMD=CCI*(TX(N)-1.)+(TX(N)+1.)
            ELSE 
               EMD=WK(N)*WK(N)*(1.-0.5*WK(N))/6.
            END IF
 
            EMP(N)=EM0*BAV+EMD*DELB
            EMM(N)=EM0*BAV-EMD*DELB
   50    CONTINUE
!
!     Accumulate partial results for clear slab.
!
         IF(INIT.EQ.-1.OR.INIT.EQ.-100)THEN
! *** First layer 
            CALL VCOPY(EMP,EMPL,NANG)
            IF(REFLG)CALL VCOPY(EMM,EMML,NANG)
            CALL VCOPY(TX,TLAY,NANG)
         ELSE 
! *** Update forward and backward emission for subsequent layers
            CALL DMXVV(TX,EMPL,WK,NANG)
            CALL sVADD(WK,EMP,EMPL,NANG)
            IF(REFLG)THEN
               CALL DMXVV(TLAY,EMM,WK,NANG)
               CALL sVADD(WK,EMML,EMML,NANG)
            END IF
! *** Update transmission
            CALL DMXVV(TLAY,TX,TLAY,NANG)
         END IF
      ELSE
!
! *** SOLAR ONLY
!
         IF(INIT.EQ.-1.OR.INIT.EQ.-100)THEN
            CALL VCOPY(WK,TLAY,NANG)
         ELSE 
            CALL sVADD(WK,TLAY,TLAY,NANG)
         END IF
!
! *** Compute transmittance for combined consecutive clear layers
!
         IF(INIT.LE.-99)THEN
            DO 80 N=1,NANG
               TLAY(N)=EXP(-TLAY(N))
  80        CONTINUE
         END IF
      END IF
      RETURN
      END subroutine clrlay
!
      SUBROUTINE MSLAY(TEXT,SSALB,TA,TB,F0,RLAY,TLAY,EMPL,EMML,&
          MAZ,VN,ALPHA,MRGUP)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the interpolation of the doubling results
!     for a cloudy layer. The routine processes NLIM spectral points.
!     Reflection, transmission and forward and backward emission are output 
!     respectively in RLAY, TLAY, EMPL and EMML. Linear-in-tau approximation 
!     is used for the thermal source and it is assumed that the Planck 
!     function varies linearly across the spectral interval spanned by the 
!     NLIM spectral points.
!
!     LAST REVISED: MARCH 1995 to implement linear spectral interpolation
!                   of the doubling results.
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL ALPHA(0:2*MXCANG)
      REAL RLAY(NANG2),TLAY(NANG2),EMPL(NANG),EMML(NANG)
      REAL EM(2*MXANG),SBM(2*MXANG)
      REAL*8 VN
      LOGICAL PLNCK,SBEAM
!
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
!
      CALL VZERO(EMPL,NANG)
      CALL VZERO(EMML,NANG)
!
      IF (PLNCK)THEN
         BA=WNPLAN(VN,TA)
         BB=WNPLAN(VN,TB)
         BAV=.5*(BB+BA)
         DELB=.5*(BB-BA)
!       Or omit linear in tau
!        bav = wnplan(vn,(ta+tb)/2.)
!        delb = 0.
      END IF
!
!     Get doubling tables for current cloud layer
!
      !CALL GENTAB(TABS,TSCAT,ALPHA,MAZ,RLAY,TLAY,EM,SBM)
      CALL GENTAB(TEXT,SSALB,ALPHA,MAZ,RLAY,TLAY,EM,SBM)
!
!     Compute forward and backward internal emission due to
!     thermal and solar sources
!     
      IF(PLNCK)THEN
!     
!     Calculate vertically averaged Planck function for the layer and
!     correction term for linear-in-Tau approximation at each spectral 
!     point
!
         DO 140 N=1,NANG
            EMPL(N)=EM(N)*BAV+EM(NANG+N)*DELB
            EMML(N)=EM(N)*BAV-EM(NANG+N)*DELB
  140    CONTINUE
      END IF
!     
      IF(SBEAM)THEN
!     
!     Add solar component
!     
         IF(MRGUP.EQ.1)THEN
            DO 230 N=1,NANG
               EMPL(N)=EMPL(N)+SBM(N)*F0
               EMML(N)=EMML(N)+SBM(NANG+N)*F0
  230       CONTINUE
         ELSE
            DO 240 N=1,NANG
               EMPL(N)=EMPL(N)+SBM(NANG+N)*F0
               EMML(N)=EMML(N)+SBM(N)*F0
  240       CONTINUE
         END IF
      END IF
!
      RETURN
      END subroutine mslay
!
      SUBROUTINE GENTAB(TEXT,SSALB,ALPHA,MAZ,R,T,EM,SBM)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     Generate look up tables of reflexion, transmission matrices
!     and source vectors by doubling calculations for the range of
!     absorption optical depths specified by NMIN and NMAX. NMIN0
!     and NMAX0 indicate the range of the calculations covered in 
!     previous calls to GENTAB so that calculations corresponding
!     to NMIN0<N<NMAX0 do not have to be redone. Table entries are 
!     in the common block /TABX/.
!
!     NOTE:
!     Internal sources due to thermal emission and solar scattering
!     are stored in vectors EM and SBM respectively. The first NANG
!     elements of SBM correspond to the reflected solar radiation,
!     while the second NANG elements correspond to the transmitted
!     solar radiation. For the thermal sources, the first NANG elements
!     contain the emission for an isothermal layer which is the same
!     in the downward and upward directions. The second set of elements 
!     are correction terms which account for the assumed linear-in-tau 
!     dependence of the PLanck function within the layer. Upward and 
!     downward thermal sources are derived in the subroutine MSLAY by 
!     respectively adding and subtracting the correction terms from the 
!     first set of elements.
!
!     REVISED:  21 MARCH 1995
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL ALPHA(0:2*MXCANG)
      REAL PFP(MXANG2),PFM(MXANG2),PFP0(MXANG),PFM0(MXANG)
      REAL R(NANG2),T(NANG2),EM(N2ANG),SBM(N2ANG)
      LOGICAL PLNCK,SBEAM
!
      COMMON /ANGL  / PTH(MXANG),GMU(MXANG),GWT(MXANG),GMU0,&
                      YLM(0:2*MXCANG-1,MXANG),YLM0(0:2*MXCANG-1)
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /CONSTN/ PI,CL2INV
!
! *** Construct phase matrix.
! *** Pre-scale matrices for diamond initialization (DI) and apply 
! *** angular weighting.
!
      DO 30 I=1,NANG
         DO 20 J=1,NANG-1
            IJ=(I-1)*NANG+J
            PFP(IJ)=0.
            PFM(IJ)=0.
            SGN=1.
            DO 10 L=MAZ,N2ANG-3
               PIJ=ALPHA(L)*YLM(L,I)*YLM(L,J)
               PFP(IJ)=PFP(IJ)+PIJ
               PFM(IJ)=PFM(IJ)+SGN*PIJ
               SGN=-SGN
  10        CONTINUE
            PFP(IJ)=-0.25*PFP(IJ)*PTH(I)*GWT(J)
            PFM(IJ)= 0.25*PFM(IJ)*PTH(I)*GWT(J)
  20     CONTINUE

! *** Observer angle
         PFP(IJ+1)=0.
         PFM(IJ+1)=0.
  30  CONTINUE
!     
! *** Construct phase matrix for solar source.
!     
      IF (SBEAM) THEN
         IF (MAZ.EQ.0) THEN
            SCAL=.25*GMU0
         ELSE
            SCAL=.5*GMU0
         END IF
!     
         DO 50 I=1,NANG
            PFP0(I)=0.
            PFM0(I)=0.
            SGN=1.
            DO 40 L=MAZ,N2ANG-3
               PI0=ALPHA(L)*YLM(L,I)*YLM0(L)
               PFP0(I)=PFP0(I)+PI0
               PFM0(I)=PFM0(I)+SGN*PI0
               SGN=-SGN
  40        CONTINUE
            PFP0(I)=SCAL*PFP0(I)*PTH(I)
            PFM0(I)=SCAL*PFM0(I)*PTH(I)
  50     CONTINUE
      ENDIF
!
! *** Use DI approximation when layer extinction optical depth
! *** is less than 2**-10. 
!
!      TEXT=TABS+TSCAT
!      SSALB=TSCAT/TEXT
! *** Compute initial layer optical depth and number of doublings
! *** necessary to reach TEXT
      EXP2=LOG(TEXT)*CL2INV
      EXP2=MAX(EXP2,-30.)
      NIT=MAX(0,INT(EXP2+4.))
      DTAU=TEXT/2.**NIT

      IF(SBEAM)ATT=EXP(-DTAU/GMU0)
! *** Initialize doubling calculations
      CALL DI(DTAU,SSALB,PFP,PFM,PFP0,PFM0,R,T,&
          EM,SBM)
      CALL DOUBL(NIT,ATT,R,T,EM,SBM)
!
      RETURN
      END subroutine gentab
!
      SUBROUTINE DI(DELTAU,SSALB,PFP,PFM,PFP0,PFM0,R,T,EM,SBM)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This routine calculates reflection, transmission and sources for
!     an infinitesimal layer using the Diamond Initialization scheme.
!     (for the treatment of the thermal emission, see note above in 
!     subroutine GENTAB)
!
!     REF: Wiscombe,1976, JQSRT, pp637-658.
!
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL PFP(NANG,NANG),PFM(NANG,NANG),PFP0(NANG),PFM0(NANG),&
           R(NANG,NANG),T(NANG,NANG),EM(NANG,2),SBM(NANG,2)
      LOGICAL PLNCK,SBEAM
!
      COMMON /ANGL  / PTH(MXANG),GMU(MXANG),GWT(MXANG),GMU0,&
                      YLM(0:2*MXCANG-1,MXANG),YLM0(0:2*MXCANG-1)
      COMMON /CONSTN/ PI,CL2INV
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / GAMMA(MXANG2),TI(MXANG2),TINV(MXANG2),&
                      RTINV(MXANG2),WK(MXANG2),WKEM(MXANG),&
                      WKSBM(2*MXANG),WK1(MXANG),WK2(MXANG)
!
      C1=DELTAU*SSALB
      C2=0.5*DELTAU
      CALL VSCAL(PFP,TI,C1,NANG2)
      CALL VSCAL(PFM,R,C1,NANG2)
      DO 10 I=1,NANG
         II=(I-1)*NANG+I
         TI(II)=TI(II)+PTH(I)*C2+1.
  10  CONTINUE
!
!     Build Gamma matrix - scale by a factor 2
!
      CALL VCOPY(TI,TINV,NANG2)
!$    CALL sMINV(TINV,NANG,NANG,WK,DET,TOL,0,1)
      CALL sMINV(TINV,NANG,DET,WK1,WK2)
      CALL MXM(TINV,NANG,R,NANG,RTINV,NANG)
      CALL MXM(R,NANG,RTINV,NANG,WK,NANG)
      CALL sVSUB(TI,WK,GAMMA,NANG2)
!$    CALL sMINV(GAMMA,NANG,NANG,WK,DET,TOL,0,1)
      CALL sMINV(GAMMA,NANG,DET,WK1,WK2)
      CALL VSCAL(GAMMA,GAMMA,2.,NANG2)
!
!     Compute reflection matrix
!
      CALL MXM(RTINV,NANG,GAMMA,NANG,R,NANG)
!
!     Compute transmission matrix
!
      CALL VCOPY(GAMMA,T,NANG2)
      DO 20 I=1,NANG
  20  T(I,I)=T(I,I)-1.
!
!     Compute thermal emission terms
!     
      IF (PLNCK)THEN
         CEM=(DELTAU-C1)*0.5
         CALL VSCAL(PTH,WKEM,CEM,NANG)
         CALL sVADD(GAMMA,R,WK,NANG2)
!        CALL MXV(WK,NANG,WKEM,NANG,EM)
         DO 40 I=1,NANG
            EM(I,1)=0.
            DO 30 J=1,NANG
               IJ=(I-1)*NANG+J
               EM(I,1)=EM(I,1)+WK(IJ)*WKEM(J) 
  30        CONTINUE
  40     CONTINUE
! *** No correction term for infinitesimal layer
         CALL VZERO(EM(1,2),NANG)
      END IF
!     
!     Compute upward and downward emission due to solar scattering
!     
      IF(SBEAM)THEN
         CSOL=.5*SSALB/PI*(1.-EXP(-DELTAU/GMU0))
         CALL VSCAL(PFM0,WKSBM,CSOL,NANG)
         CALL VSCAL(PFP0,WKSBM(NANG+1),CSOL,NANG)
         DO 60 I=1,NANG
            SBM(I,1)=0.
            SBM(I,2)=0.
            DO 50 J=1,NANG
               IJ=(I-1)*NANG+J
               SBM(I,1)=SBM(I,1)+GAMMA(IJ)*WKSBM(J)+R(J,I)*WKSBM(NANG+J)
               SBM(I,2)=SBM(I,2)+GAMMA(IJ)*WKSBM(NANG+J)+R(J,I)*WKSBM(J)
  50        CONTINUE                                                    
  60     CONTINUE                                                       
      END IF
      RETURN
      END subroutine di
!
      SUBROUTINE DOUBL(NIT,ATT,R,T,EM,SBM)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the doubling calculations for a homogeneous 
!     scattering layer. The number of calculations is specified by NIT.
!     It is assumed that the reflection and transmission matrices and the
!     source vectors have already been initialized.
!
!     REVISED:  21 MARCH 1995
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL R(NANG2),T(NANG2),EM(N2ANG),SBM(N2ANG)
      LOGICAL PLNCK,SBEAM
!
      COMMON /LPAR  / PLNCK,SBEAM
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / RM(MXANG2),WK1(MXANG2),WK2(MXANG2),WKT(MXANG2),&
                      WKEM(2*MXANG) 
!
      DO 20 N=1,NIT
!
!     Compute multiple reflection term
!
         CALL MXM(R,NANG,R,NANG,RM,NANG)
         DO 10 I=1,NANG2
            RM(I)=-RM(I)
  10     CONTINUE
         CALL MADDI(RM,NANG)
!$       CALL sMINV(RM,NANG,NANG,WK1,DET,TOL,0,1)
         CALL sMINV(RM,NANG,DET,WK1,WK2)
! 
!     Compute new transmission matrix
!
         CALL MXM(RM,NANG,T,NANG,WK1,NANG)
         CALL VCOPY(T,WKT,NANG2)
         CALL MXM(WKT,NANG,WK1,NANG,T,NANG)
!
!     Double thermal and solar sources
!     
         CALL MXM(R,NANG,WK1,NANG,WK2,NANG)
         IF(PLNCK)CALL EMIS(EM,WK1,WK2,WKEM)
         IF(SBEAM)CALL SUNB(SBM,WK1,WK2,WKEM,ATT)
!     
!     Compute new reflection matrix
!
         CALL MXM(WKT,NANG,WK2,NANG,WK1,NANG)
         CALL sVADD(R,WK1,R,NANG2)
  20  CONTINUE
      RETURN
      END subroutine doubl
!
      SUBROUTINE EMIS(EM,TM,TMR,WKEM)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the doubling of the thermal source
!     (see note in subroutine GENTAB)
!
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL EM(NANG,2),WKEM(NANG,2),TM(NANG,NANG),TMR(NANG,NANG)
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
!
      CALL VCOPY(EM,WKEM,N2ANG)
      DO 20 I=1,NANG
         EM(I,2)=EM(I,2)+EM(I,1)
         DO 10 J=1,NANG
            EM(I,1)=EM(I,1)+(TM(J,I)+TMR(J,I))*WKEM(J,1)
            EM(I,2)=EM(I,2)+(TM(J,I)-TMR(J,I))*(WKEM(J,2)-WKEM(J,1))
  10     CONTINUE
  20  CONTINUE
      CALL VSCAL(EM(1,2),EM(1,2),.5,NANG)
      RETURN
      END subroutine emis
!
      SUBROUTINE SUNB(SBM,TM,TMR,WKEM,ATT)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the doubling of the solar source
!     (see note in subroutine GENTAB)
!
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL SBM(NANG,2),WKEM(NANG,2),TM(NANG,NANG),TMR(NANG,NANG)
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
!
      CALL VCOPY(SBM,WKEM,N2ANG)
      CALL VSCAL(SBM(1,2),SBM(1,2),ATT,NANG)
      CALL VSCAL(WKEM,WKEM,ATT,NANG)
      DO 20 I=1,NANG
         DO 10 J=1,NANG
            SBM(I,1)=SBM(I,1)+TM(J,I)*WKEM(J,1)+TMR(J,I)*WKEM(J,2)
            SBM(I,2)=SBM(I,2)+TM(J,I)*WKEM(J,2)+TMR(J,I)*WKEM(J,1)
  10     CONTINUE
  20  CONTINUE
      ATT=ATT*ATT
      RETURN
      END subroutine sunb
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     OSS (Optimal Spectral Sampling)
!     $Name:  $ 
!     $Id: oss_addbl.f90,v 1.4 2006/12/15 21:04:57 raschbre Exp $ 
!     Copyright AER, Inc., 2002, 2003. All rights Reserved.
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE ADDCLR(R,T,EM,REFLG,PLNCK)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine does the adding of a clear layer - Transmission
!     and internal emission for the layer must be contained in common
!     block /LAYDAT/. Results are accumulated in arrays R, T and EM.
!     Transmittance is updated only if TXON = .TRUE., and reflection
!     matrix only if non-zero initially (REFLG=.FALSE.).
!
!     REVISED:  08 MARCH 1993  
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     
      REAL R(*),T(*),EM(*)
      LOGICAL REFLG,PLNCK
!     
      COMMON /CONSTN/ PI,CL2INV
      COMMON /LAYDAT/ RLAY(MXANG2),TLAY(MXANG2),&
                      EMPL(MXANG),EMML(MXANG)
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / WKEM(MXANG),TR(MXANG2),&
                      WKT(MXANG2)
!
!     Update emission
!
      CALL DMXVV(TLAY,EM,WKEM,NANG)
!
      IF(PLNCK)THEN
         CALL SVADD(WKEM,EMPL,EM,NANG)
      ELSE
         CALL VCOPY(WKEM,EM,NANG)
      END IF
!
      IF(REFLG)THEN
!
         CALL DMXMV(TLAY,R,TR,NANG)
!
         IF(PLNCK)THEN
            CALL MXVV(TR,EMML,WKEM,NANG)
            CALL SVADD(WKEM,EM,EM,NANG)
         END IF
!
!     Update reflection matrix
!
         CALL MXDMV(TR,TLAY,R,NANG)
      END IF

      RETURN
      END subroutine addclr
!
      SUBROUTINE ADDMS(R,T,EM,REFLG)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     General layer adding routine - Reflection, transmission and 
!     internal emission for the layer must be contained in common
!     block /LAYDAT/. Results are accumulated in arrays R, T and EM.
!     Transmittance is updated only if TXON = .TRUE., and reflection
!     matrix only if non-zero initially (REFLG=.FALSE.).
!
!     REVISED:  08 MARCH 1993  
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      LOGICAL REFLG
      REAL R(*),T(*),EM(*)
!
      COMMON /LAYDAT/ RL(MXANG2),TL(MXANG2),&
                      EMPL(MXANG),EMML(MXANG)
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / RM(MXANG2),WK1(MXANG2),WK2(MXANG2),WKEM1(MXANG), WKEM2(MXANG)
!
!      EQUIVALENCE (WK0,RM) 
!
      IF(REFLG)THEN
!
!     Calculate multiple reflection term (result in RM)
!     
         CALL MLREF(R,RL)
!     
!     Update emission
!
         CALL MXMV(TL,RM,WK1,NANG)
         CALL MXMV(WK1,R,WK2,NANG)
         CALL MXVV(WK1,EM,WKEM1,NANG)
         CALL MXVV(WK2,EMML,WKEM2,NANG)
         CALL SVADD(WKEM1,WKEM2,EM,NANG)
         CALL SVADD(EM,EMPL,EM,NANG)
!     
!     Update reflection matrix
!
!       -- old call:
!         CALL MXMV(WK2,TL,WK0,NANG)
!         CALL SVADD(WK0,RL,R,NANG2)
!       -- new call, removing equivalence:
         CALL MXMV(WK2,TL,RM,NANG)
         CALL SVADD(RM,RL,R,NANG2)
!
      ELSE
!
         CALL MXVV(TL,EM,WKEM1,NANG)
         CALL SVADD(EMPL,WKEM1,EM,NANG)
!
         CALL VCOPY(RL,R,NANG2)
!
         REFLG=.TRUE.
!
      END IF
      RETURN
      END subroutine addms
!
      SUBROUTINE MRGRAD(RF,EMF,RB,EMB,REFLF,REFLB,RAD)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     This subroutine calculates the observed radiances by combining 
!     reflection and emission of the atmosphere below and above the 
!     observer level.
!     
!     Flags REFLF/REFLB indicate whether the portion of the atmosphere
!     in front of/behind the observer is reflecting
!
!     REVISED:  11 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL RF(*),EMF(*),RB(*),EMB(*),RAD(*)
      LOGICAL REFLF,REFLB
!
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / RM(MXANG2),RADTMP(MXANG)
!
      IF (.NOT.REFLF)THEN
         CALL VCOPY(EMF,RAD,NANG)
      ELSE IF (.NOT.REFLB)THEN
         CALL MXVV(RF,EMB,RAD,NANG)
         CALL sVADD(EMF,RAD,RAD,NANG)
      ELSE
         CALL MLREF(RF,RB)
         CALL MXVV(RF,EMB,RADTMP,NANG)
         CALL sVADD(EMF,RADTMP,RADTMP,NANG)
         CALL MXVV(RM,RADTMP,RAD,NANG)
      END IF
      RETURN
      end subroutine mrgrad
!
      SUBROUTINE MLREF(RA,RB)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     Calculates multiple reflection term by truncated series.
!
!     REVISED:  11 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL RA(*),RB(*)
!
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      COMMON /LOCWK / RM(MXANG2),RAB(MXANG2), Q1(MXANG2),Q2(MXANG2)
!
      CALL MXMV(RA,RB,RAB,NANG)
!
      CALL MXMV(RAB,RAB,Q1,NANG)
      CALL SVADD(RAB,Q1,RM,NANG2)
!
      CALL MXMV(RAB,Q1,Q2,NANG)
      CALL SVADD(RM,Q2,RM,NANG2)
!
      CALL MXMV(RAB,Q2,Q1,NANG)
      CALL SVADD(RM,Q1,RM,NANG2)
!
      CALL MXMV(RAB,Q1,Q2,NANG)
      CALL HOCOR(RM,Q1,Q2)
!
      CALL MADDI(RM,NANG)
!
      RETURN
      end subroutine mlref
!
      SUBROUTINE HOCOR(RM,Q1,Q2)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!     Correction for missing high order terms in MLREF
!
!     REVISED:  13 APRIL 1994
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL RM(*),Q1(*),Q2(*)
!
      COMMON /PNLPAR/ NANG,N2ANG,NANG2
      DATA EPSIL/1.E-16/
!
      DO 20 L=1,NANG2
         IF(MOD(L,NANG).NE.0)THEN
            RM(L)=RM(L)+Q2(L)*Q1(L)/((Q1(L)-Q2(L))+EPSIL)
         END IF
  20  CONTINUE
      RETURN
      end subroutine hocor
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     OSS (Optimal Spectral Sampling)
!     $Name:  $ 
!     $Id: oss_addbl.f90,v 1.4 2006/12/15 21:04:57 raschbre Exp $ 
!     Copyright AER, Inc., 2002, 2003. All rights Reserved.
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     VECTOR OPERATIONS     
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE VZERO(A,N)
! **********************************************************************
! *   Set vector A to zero.                                            *
! **********************************************************************
      DIMENSION A(N)
      DO 10 I=1,N
   10 A(I)=0.0
      RETURN
      end subroutine vzero
!
      SUBROUTINE VCOPY(A,B,N)
! **********************************************************************
! *   Copy vector A to vector B.                                       *
! **********************************************************************
      DIMENSION A(N),B(N)
      DO 10 I=1,N
   10 B(I)=A(I)
      RETURN
      end subroutine vcopy
!
      SUBROUTINE sVADD(A,B,C,N)
! **********************************************************************
! *   Add vectors A and B.  Result is vector C                         *
! **********************************************************************
      DIMENSION A(N),B(N),C(N)
      DO 10 I=1,N
   10 C(I)=A(I)+B(I)
      RETURN
      end subroutine sVADD
!
      SUBROUTINE VSCAL(A,B,S,N)
! **********************************************************************
! *   Multiply vector A by scalar S. Result is vector B.               *
! **********************************************************************
      DIMENSION A(N),B(N)
      DO 10 I=1,N
   10 B(I)=A(I)*S
      RETURN
      end subroutine vscal
!
      SUBROUTINE sVSUB(A,B,C,N)
! **********************************************************************
! *   Subtract vector B from vector A.  Result is vector C.            *
! **********************************************************************
      DIMENSION A(N),B(N),C(N)
      DO 10 I=1,N
   10 C(I)=A(I)-B(I)
      RETURN
      END subroutine svsub
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     MATRIX AND MATRIX/VECTOR OPERATIONS 
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE MADDI(A,ISIZ)
! **********************************************************************
! *   Add the identity matrix to matrix A.                             *
! **********************************************************************
      DIMENSION A(ISIZ*ISIZ)
      DO 10 I=1,ISIZ*ISIZ,ISIZ+1
   10 A(I)=A(I)+1.
      RETURN
      END subroutine maddi
!
      SUBROUTINE MXM(A,NAR,B,NAC,C,NBC)
! **********************************************************************
! *   CALCULATES PRODUCT OF TWO MATRICES--C=BA AND ASSUMES THE         *
! *   SKIP DISTANCE BETWEEN ELEMENTS TO BE 1.                          *
! *                                                                    *
! *   REVISED:  11-11-93 (D. PORTMAN)                                  *
! *                                                                    *
! *   ARGUMENTS:                                                       *
! *                                                                    *
! *   A       (input) Second matrix of product.                        *
! *   NAR     (input/output) Number of rows of matrices A and C.       *
! *   B       (input) First matrix of product.                         *
! *   NAC     (input) Number of columns of matrix A and the number of  *
! *           rows of matrix B.                                        *
! *   C       (output) Result matrix.                                  *
! *   NBC     (input/output) Number of columns of matrices B and C.    *
! *                                                                    *
! **********************************************************************
      REAL A(NAR,NAC), B(NAC,NBC), C(NAR,NBC)
!
      DO 30 I = 1, NAR
          DO 20 J = 1, NBC
              SUM = 0.
              DO 10 K = 1, NAC 
                  SUM = SUM + (A(I,K)*B(K,J))
   10         CONTINUE
              C(I,J) = SUM
   20     CONTINUE
   30 CONTINUE
      RETURN
      end subroutine mxm
!
      SUBROUTINE sMINV(A,N,D,L,M)
! **********************************************************************
! ROUTINE          DMINV   SUBROUTINE  *****  [EYRE.RETCOF]
!
! PURPOSE          TO INVERT A DOUBLE PRECISION MATRIX
!
! VERSION          3.00,150385,J.R.EYRE
!
! DESCRIPTION      ROUTINE TO INVERT A DOUBLE PRECISION MATRIX.
!                  METHOD: THE STANDARD GAUSS-JORDAN METHOD IS USED.
!                  THE DETERMINANT IS ALSO CALCULATED. A DETERMINANT
!                  OF ZERO INDICATES THAT THE MATRIX IS SINGULAR.
!
! ARGUMENTS        A      R*8 INPUT MATRIX, DESTROYED IN COMPUTATION AND
!                             REPLACED BY RESULTANT INVERSE.
!                  N      I*4 ORDER OF MATRIX A
!                  D      R*8 RESULTANT DETERMINANT
!                  L      R*8 WORK VECTOR OF LENGTH N
!                  M      R*8 WORK VECTOR OF LENGTH N
!
!     ..................................................................
!
      DIMENSION A(N*N),L(N),M(N)
!
!     ..................................................................
!
!        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
!        C$ IN COLUMNS 1-2 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
!        STATEMENT WHICH FOLLOWS.
!
!$      REAL*8 A,D,BIGA,HOLD
        REAL A,D,BIGA,HOLD,L,M     
!
!        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
!        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
!        ROUTINE.
!
!        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
!        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  ABS IN STATEMENT
!        10 MUST BE CHANGED TO DABS.
!
! **********************************************************************
!
!        SEARCH FOR LARGEST ELEMENT
!
      D=1.0
      NK=-N
      DO 80 K=1,N
         NK=NK+N
         L(K)=K
         M(K)=K
         KK=NK+K
         BIGA=A(KK)
         DO 20 J=K,N
             IZ=N*(J-1)
             DO 20 I=K,N
                IJ=IZ+I
!$ 10           IF(DABS(BIGA)-DABS(A(IJ))) 15,20,20
   10           IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20    
   15           BIGA=A(IJ)
                L(K)=I
                M(K)=J
   20    CONTINUE
!
!        INTERCHANGE ROWS
!
         J=L(K)
         IF(J-K) 35,35,25
   25    KI=K-N
         DO 30 I=1,N
            KI=KI+N
            HOLD=-A(KI)
            JI=KI-K+J
            A(KI)=A(JI)
            A(JI) =HOLD
   30    CONTINUE 
!
!        INTERCHANGE COLUMNS
!
   35    I=M(K)
         IF(I-K) 45,45,38
   38    JP=N*(I-1)
         DO 40 J=1,N
            JK=NK+J
            JI=JP+J
            HOLD=-A(JK)
            A(JK)=A(JI)
            A(JI) =HOLD
   40    CONTINUE
!
!        DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
!        CONTAINED IN BIGA)
!
   45    IF(BIGA) 48,46,48
   46    D=0.0
         RETURN
   48    DO 55 I=1,N
            IF(I-K) 50,55,50
   50       IK=NK+I
            A(IK)=A(IK)/(-BIGA)
   55    CONTINUE
!
!        REDUCE MATRIX
!
         DO 65 I=1,N
            IK=NK+I
            HOLD=A(IK)
            IJ=I-N
            DO 65 J=1,N
               IJ=IJ+N
               IF(I-K) 60,65,60
   60          IF(J-K) 62,65,62
   62          KJ=IJ-I+K
               A(IJ)=HOLD*A(KJ)+A(IJ)
   65    CONTINUE
!
!        DIVIDE ROW BY PIVOT
!
         KJ=K-N
         DO 75 J=1,N
            KJ=KJ+N
            IF(J-K) 70,75,70
   70       A(KJ)=A(KJ)/BIGA
   75    CONTINUE
!
!        PRODUCT OF PIVOTS
!
         D=D*BIGA
!
!        REPLACE PIVOT BY RECIPROCAL
!
         A(KK)=1.0/BIGA
   80 CONTINUE
!
!        FINAL ROW AND COLUMN INTERCHANGE
!
      K=N
  100 K=(K-1)
      IF(K) 150,150,105
  105 I=L(K)
      IF(I-K) 120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
         JK=JQ+J
         HOLD=A(JK)
         JI=JR+J
         A(JK)=-A(JI)
         A(JI) =HOLD
  110 CONTINUE
  120 J=M(K)
      IF(J-K) 100,100,125
  125 KI=K-N
      DO 130 I=1,N
         KI=KI+N
         HOLD=A(KI)
         JI=KI-K+J
         A(KI)=-A(JI)
         A(JI) =HOLD
  130 CONTINUE
      GO TO 100
  150 RETURN
      END subroutine sminv
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     VECTORIZED CHAIN MATRIX AND MATRIX/VECTOR OPERATIONS    
!     REVISED:  04 OCTOBER 1993
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE MXMV(A,B,C,ISIZ)
      DIMENSION A(ISIZ,ISIZ),B(ISIZ,ISIZ),C(ISIZ,ISIZ)
      DO 50 I=1,ISIZ
         DO 40 J=1,ISIZ
            C(J,I)=0.
            DO 30 K=1,ISIZ
               C(J,I)=C(J,I)+A(K,I)*B(J,K)
  30        CONTINUE
  40     CONTINUE
  50  CONTINUE
      RETURN
      END subroutine mxmv
!
      SUBROUTINE DMXMV(A,B,C,ISIZ)
      DIMENSION A(ISIZ),B(ISIZ,ISIZ),C(ISIZ,ISIZ)
      DO 30 I=1,ISIZ
         DO 20 J=1,ISIZ
            C(J,I)=A(I)*B(J,I)
  20     CONTINUE
  30  CONTINUE
      RETURN
      END subroutine dmxmv
!
      SUBROUTINE MXDMV(A,B,C,ISIZ)
      DIMENSION A(ISIZ,ISIZ),B(ISIZ),C(ISIZ,ISIZ)
      DO 30 I=1,ISIZ
         DO 20 J=1,ISIZ
            C(J,I)=A(J,I)*B(J)
  20     CONTINUE
  30  CONTINUE
      RETURN
      END subroutine mxdmv
!
      SUBROUTINE MXVV(A,B,C,ISIZ)
      DIMENSION A(ISIZ,ISIZ),B(ISIZ),C(ISIZ)
      DO 40 I=1,ISIZ
         C(I)=0.
         DO 30 J=1,ISIZ
            C(I)=C(I)+A(J,I)*B(J)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END subroutine mxvv
!
      SUBROUTINE DMXVV(A,B,C,ISIZ)
      DIMENSION A(ISIZ),B(ISIZ),C(ISIZ)
      DO 20 I=1,ISIZ
         C(I)=A(I)*B(I)
   20 CONTINUE
      RETURN
      END subroutine dmxvv


      SUBROUTINE  LEPOLY_OSS( NMU, M, MAXMU, TWONM1, MU, YLM )
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!       ORIGIN:     DISORT (Stamnes etal.) V1.1
!
!       COMPUTES THE NORMALIZED ASSOCIATED LEGENDRE POLYNOMIAL,
!       DEFINED IN TERMS OF THE ASSOCIATED LEGENDRE POLYNOMIAL
!       PLM = P-SUB-L-SUPER-M AS
!
!             YLM(MU) = SQRT( (L-M)!/(L+M)! ) * PLM(MU)
!
!       FOR FIXED ORDER -M- AND ALL DEGREES FROM L = M TO TWONM1.
!       WHEN M.GT.0, ASSUMES THAT Y-SUB(M-1)-SUPER(M-1) IS AVAILABLE
!       FROM A PRIOR CALL TO THE ROUTINE.
!
!       REFERENCE: Dave, J.V. and B.H. Armstrong, Computations of
!                  High-Order Associated Legendre Polynomials,
!                  J. Quant. Spectrosc. Radiat. Transfer 10,
!                  557-562, 1970.  (hereafter D/A)
!
!       METHOD: Varying degree recurrence relationship.
!
!       NOTE 1: The D/A formulas are transformed by
!               setting  M = n-1; L = k-1.
!       NOTE 2: Assumes that routine is called first with  M = 0,
!               then with  M = 1, etc. up to  M = TWONM1.
!       NOTE 3: Loops are written in such a way as to vectorize.
!
!  I N P U T     V A R I A B L E S:
!
!       NMU    :  NUMBER OF ARGUMENTS OF -YLM-
!       M      :  ORDER OF -YLM-
!       MAXMU  :  FIRST DIMENSION OF -YLM-
!       TWONM1 :  MAX DEGREE OF -YLM-
!       MU(I)  :  I = 1 TO NMU, ARGUMENTS OF -YLM-
!       IF M.GT.0, YLM(M-1,I) FOR I = 1 TO NMU IS REQUIRED
!
!  O U T P U T     V A R I A B L E:
!
!       YLM(L,I) :  L = M TO TWONM1, NORMALIZED ASSOCIATED LEGENDRE
!                   POLYNOMIALS EVALUATED AT ARGUMENT -MU(I)-
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL  MU(*), YLM( 0:MAXMU,* )
      INTEGER M, NMU, TWONM1
      PARAMETER ( MAXSQT = 1000 )
      REAL SQT( MAXSQT )
      LOGICAL PASS1
      SAVE SQT, PASS1
      DATA PASS1 / .TRUE. /
!
      IF ( PASS1 )  THEN
         PASS1 = .FALSE.
         DO 10 NS = 1, MAXSQT
            SQT( NS ) = SQRT( FLOAT(NS) )
   10    CONTINUE
      ENDIF
!
      IF ( 2*TWONM1 .GT. MAXSQT ) STOP  'LEPOLY--NEED TO INCREASE PARAM MAXSQT'
!
      IF ( M .EQ. 0 )  THEN
!
!     UPWARD RECURRENCE FOR ORDINARY LEGENDRE POLYNOMIALS
!
         DO  20  I = 1, NMU
            YLM( 0,I ) = 1.
            YLM( 1,I ) = MU( I )
   20    CONTINUE
!
         DO  40  L = 2, TWONM1
            DO  30  I = 1, NMU
               YLM( L,I ) = ( ( 2*L-1 ) * MU(I) * YLM( L-1,I ) &
                            - ( L-1 ) * YLM( L-2,I ) ) / L
   30       CONTINUE
   40   CONTINUE
      ELSE
         DO  50  I = 1, NMU
!
!     Y-SUB-M-SUPER-M; DERIVED FROM D/A EQS. (11,12)
!
            YLM( M,I) = - SQT( 2*M-1 ) / SQT( 2*M ) &
                        * SQRT( 1. - MU(I)**2 ) * YLM( M-1,I )
!
!     Y-SUB-(M+1)-SUPER-M; DERIVED FROM D/A EQS. (13,14) USING EQS. (11,12)
!
            YLM( M+1,I ) = SQT( 2*M+1 ) * MU(I) * YLM( M,I )
   50    CONTINUE
!
!     UPWARD RECURRENCE; D/A EQ. (10)
!
         DO  70  L = M+2, TWONM1
            TMP1 = SQT( L-M ) * SQT( L+M )
            TMP2 = SQT( L-M-1 ) * SQT( L+M-1 )
            DO  60  I = 1, NMU
               YLM( L,I ) = ( ( 2*L-1 ) * MU(I) * YLM( L-1,I ) &
                              - TMP2 * YLM( L-2,I ) ) / TMP1
   60       CONTINUE
   70   CONTINUE
      END IF
      RETURN
      END subroutine lepoly_oss




!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     OSS (Optimal Spectral Sampling)
!     $Name:  $ 
!     $Id: oss_addbl.f90,v 1.4 2006/12/15 21:04:57 raschbre Exp $ 
!     Copyright AER, Inc., 2002, 2003. All rights Reserved.
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      function wdraddt(vn,t)
      IMPLICIT REAL*8 (A-H,O-Z)
      real t,wdraddt
      parameter(h=6.626176d-27,c=2.997925d+10,b=1.380662d-16)
      parameter(c1=2.d0*h*c*c)
      parameter(c2=h*c/b)
!
      t_1=1/t
!
      ct1v1=c1*vn**3
      ct2v1=c2*vn*t_1
      ct3v1=exp(ct2v1)
!
      bplan=ct1v1/(ct3v1-1.)
      wdraddt=bplan/(ct3v1-1.)*ct2v1*ct3v1*t_1
!
      return
      end function wdraddt

      FUNCTION WNBRIT(VN,RAD)
! $ RADIANCE TO BRIGHTNESS TEMPERATURE (HMW)
! * 'NEW' PLANCK'S CONSTANT, VELOCITY OF LIGHT, BOLTZMANN'S CONSTANT
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL RAD,WNBRIT
      PARAMETER (H = 6.626176D-27, C = 2.997925D+10, B = 1.380662D-16)
      PARAMETER (C1 = 2.D0*H*C*C)
      PARAMETER (C2 = H*C/B)
      TNB(X,Y,Z)=Y/DLOG(X/Z+1.D0)
!
      F1=C1*VN**3
      F2=C2*VN
      R=RAD
      TBB=TNB(F1,F2,R)
      WNBRIT=TBB
      RETURN
      END function wnbrit

      FUNCTION DBDT(VN,T)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL T,DBDT
      PARAMETER(H=6.626176D-27,C=2.997925D+10,B=1.380662D-16)
      PARAMETER(C1=2.D0*H*C*C)
      PARAMETER(C2=H*C/B)
!
      T_1=1/T
!
      CT1V1=C1*VN**3
      CT2V1=C2*VN*T_1
      CT3V1=EXP(CT2V1)
!
      BPLAN=CT1V1/(CT3V1-1.)
      DBDT=BPLAN/(CT3V1-1.)*CT2V1*CT3V1*T_1
!
      RETURN
      END function dbdt
                                                                   
      FUNCTION WNMRAD(VN,TAU,TEM,TS,LS)
! $ MONOCHROMATIC RADIANCE CALCULATION (HMW)
! Units are mw*cm/m2/sr
      REAL*8 VN
      DIMENSION TAU(*),TEM(*)
      TAU1=TAU(1)
      T1=TEM(1)
      B1=WNPLAN(VN,T1)
      RAD=0.
      DO 120 I=2,LS
      TAU2=TAU(I)
      T2=TEM(I)
      B2=WNPLAN(VN,T2)
      RAD=RAD+.5*(B1+B2)*(TAU1-TAU2)
      TAU1=TAU2
  120 B1=B2
      BSTS=WNPLAN(VN,TS)*TAU1
      RAD=RAD+BSTS
      WNMRAD=RAD
      RETURN
      END function wnmrad

      FUNCTION WNPLAN(VN,TEM)
! $ TEMPERATURE TO PLANCK RADIANCE (HMW)
! $ 'NEW' PLANCK'S CONSTANT, VELOCITY OF LIGHT, BOLTZMANN'S CONSTANT
! Units are mw*cm/m2/sr
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL TEM,WNPLAN
      PARAMETER (H = 6.626176D-27, C = 2.997925D+10, B = 1.380662D-16)
! Units of H are mw, units of c are cm/s
      PARAMETER (C1 = 2.D0*H*C*C)
      PARAMETER (C2 = H*C/B)
      BNT(X,Y,Z)=X/(DEXP(Y/Z)-1.D0)
!
      !---To avoid dividing by zero
      IF (TEM.EQ.0.) THEN
         WNPLAN=0.
         RETURN
      ENDIF
      F1=C1*VN**3
      F2=C2*VN
      T=TEM
      RAD=BNT(F1,F2,T)
      WNPLAN=RAD
      !---test
      !WNPLAN=TEM
      !------
      RETURN
      END function wnplan

      SUBROUTINE QGAUSN( M, GMU, GWT )

!       Compute weights and abscissae for ordinary Gaussian quadrature
!       on the interval (0,1);  that is, such that

!           sum(i=1 to M) ( GWT(i) f(GMU(i)) )

!       is a good approximation to

!           integral(0 to 1) ( f(x) dx )

!   INPUT :    M       order of quadrature rule

!   OUTPUT :  GMU(I)   array of abscissae (I = 1 TO M)
!             GWT(I)   array of weights (I = 1 TO M)

!   REFERENCE:  Davis, P.J. and P. Rabinowitz, Methods of Numerical
!                   Integration, Academic Press, New York, pp. 87, 1975

!   METHOD:  Compute the abscissae as roots of the Legendre
!            polynomial P-sub-M using a cubically convergent
!            refinement of Newton's method.  Compute the
!            weights from EQ. 2.7.3.8 of Davis/Rabinowitz.  Note
!            that Newton's method can very easily diverge; only a
!            very good initial guess can guarantee convergence.
!            The initial guess used here has never led to divergence
!            even for M up to 1000.

!   ACCURACY:  relative error no better than TOL or computer
!              precision (machine epsilon), whichever is larger

!   INTERNAL VARIABLES:

!    ITER      : number of Newton Method iterations
!    MAXIT     : maximum allowed iterations of Newton Method
!    PM2,PM1,P : 3 successive Legendre polynomials
!    PPR       : derivative of Legendre polynomial
!    P2PRI     : 2nd derivative of Legendre polynomial
!    TOL       : convergence criterion for Legendre poly root iteration
!    X,XI      : successive iterates in cubically-convergent version
!                of Newtons Method (seeking roots of Legendre poly.)

!   Called by- SETDIS, SURFAC
!   Calls- D1MACH, ERRMSG
! +-------------------------------------------------------------------+

!     .. Scalar Arguments ..

      INTEGER   M
!     ..
!     .. Array Arguments ..

      REAL      GMU( M ), GWT( M )
!     ..
!     .. Local Scalars ..

      INTEGER   ITER, K, LIM, MAXIT, NN, NP1
      REAL      CONA, PI, T
      DOUBLE PRECISION EN, NNP1, ONE, P, P2PRI, PM1, PM2, PPR, PROD,&
                      TMP, TOL, TWO, X, XI
!     ..
!     .. Intrinsic Functions ..

      INTRINSIC ABS, ASIN, COS, FLOAT, MOD, TAN
!     ..
      SAVE      PI, TOL

      DATA      PI / 0.0 / , MAXIT / 1000 / , ONE / 1.D0 / , TWO / 2.D0 /


      IF( PI.EQ.0.0 ) THEN

         PI   = 2.*ASIN( 1.0 )
         TOL  = 10.*D1MACH( 4 )

      END IF


      IF( M.LT.1 ) CALL ERRMSG( 'QGAUSN--Bad value of M',.True.)

      IF( M.EQ.1 ) THEN

         GMU( 1 ) = 0.5
         GWT( 1 ) = 1.0
         RETURN

      END IF

      EN   = M
      NP1  = M + 1
      NNP1 = M*NP1
      CONA = FLOAT( M - 1 ) / ( 8*M**3 )

      LIM  = M / 2

      DO 30 K = 1, LIM
!                                        ** Initial guess for k-th root
!                                        ** of Legendre polynomial, from
!                                        ** Davis/Rabinowitz (2.7.3.3a)
         T  = ( 4*K - 1 )*PI / ( 4*M + 2 )
         X  = COS( T + CONA / TAN( T ) )
         ITER = 0
!                                        ** Upward recurrence for
!                                        ** Legendre polynomials
   10    CONTINUE
         ITER   = ITER + 1
         PM2    = ONE
         PM1    = X

         DO 20 NN = 2, M
            P    = ( ( 2*NN - 1 )*X*PM1 - ( NN - 1 )*PM2 ) / NN
            PM2  = PM1
            PM1  = P
   20    CONTINUE
!                                              ** Newton Method
         TMP    = ONE / ( ONE - X**2 )
         PPR    = EN*( PM2 - X*P )*TMP
         P2PRI  = ( TWO*X*PPR - NNP1*P )*TMP
         XI     = X - ( P / PPR )*( ONE +( P / PPR )*P2PRI / ( TWO*PPR ) )

!                                              ** Check for convergence
         IF( ABS( XI - X ).GT.TOL ) THEN

            IF( ITER.GT.MAXIT ) &
                CALL ERRMSG( 'QGAUSN--max iteration count',.True.)

            X  = XI
            GO TO  10

         END IF
!                             ** Iteration finished--calculate weights,
!                             ** abscissae for (-1,1)
         GMU( K ) = -X
         GWT( K ) = TWO / ( TMP*( EN*PM2 )**2 )
         GMU( NP1 - K ) = -GMU( K )
         GWT( NP1 - K ) = GWT( K )
   30 CONTINUE
!                                    ** Set middle abscissa and weight
!                                    ** for rules of odd order
      IF( MOD( M,2 ).NE.0 ) THEN

         GMU( LIM + 1 ) = 0.0
         PROD   = ONE

         DO 40 K = 3, M, 2
            PROD   = PROD * K / ( K - 1 )
   40    CONTINUE

         GWT( LIM + 1 ) = TWO / PROD**2
      END IF

!                                        ** Convert from (-1,1) to (0,1)
      DO 50 K = 1, M
         GMU( K ) = 0.5*GMU( K ) + 0.5
         GWT( K ) = 0.5*GWT( K )
   50 CONTINUE


      RETURN
      END subroutine qgausn

      SUBROUTINE  ErrMsg( MESSAG, FATAL )

!        Print out a warning or error message;  abort if error

      LOGICAL       FATAL, MsgLim
      CHARACTER*(*) MESSAG
      INTEGER       MaxMsg, NumMsg
      SAVE          MaxMsg, NumMsg, MsgLim
      DATA NumMsg / 0 /,  MaxMsg / 100 /,  MsgLim / .FALSE. /


      IF ( FATAL )  THEN
         WRITE ( *, '(/,2A,/)' )  ' ******* ERROR >>>>>>  ', MESSAG
         STOP
      END IF

      NumMsg = NumMsg + 1
      IF( MsgLim )  RETURN

      IF ( NumMsg.LE.MaxMsg )  THEN
         WRITE ( *, '(/,2A,/)' )  ' ******* WARNING >>>>>>  ', MESSAG
      ELSE
         WRITE ( *,99 )
         MsgLim = .True.
      ENDIF

      RETURN

   99 FORMAT( //,' >>>>>>  TOO MANY WARNING MESSAGES --  ', &
         'They will no longer be printed  <<<<<<<', // )
 
      END subroutine errmsg

      DOUBLE PRECISION FUNCTION D1MACH(I)

!  Double-precision machine constants (see R1MACH for documentation).

!  By default, returns values appropriate for a computer with IEEE 
!  arithmetic.  This is an abbreviated version of a routine widely
!  used for 20+ years by numerical analysts.  Most of the values in
!  the original version pertain to computers which went to computer
!  heaven years ago and are of little if any interest.
! 
!  If the values herein do not work for any reason, just look in
!  your Fortran manual for the correct values (usually in the part
!  discussing representations of numbers) and insert them. The exact
!  values are not that important; they can be a factor of 2-3 off
!  without causing any harm.

!  Only I = 1,2,4 is actually used by DISORT. 

!  This routine is superseded in Fortran-90 by the intrinsic numeric 
!  inquiry functions HUGE(1.D0), TINY(1.D0), and EPSILON(1.D0).

!  The original version can be found on NetLib (search by name):
!      http://www.netlib.org/
! ====================================================================

      INTEGER   I
      !EXTERNAL  ERRMSG

      IF( I.EQ.1 )  THEN
         D1MACH = 2.3D-308
!        D1MACH = TINY(1.D0)
      ELSE IF( I.EQ.2 )  THEN  
         D1MACH = 1.7D+308
!        D1MACH = HUGE(1.D0)
      ELSE IF( I.EQ.4 )  THEN  
         D1MACH = 2.3D-16
!        D1MACH = EPSILON(1.D0)
      ELSE
         CALL ERRMSG( 'D1MACH--argument incorrect', .TRUE.)
      END IF

      RETURN
      END function d1mach

  end module oss_addbl