text:
! ###########################################################
! #                                                         #
! #                    THE LIDORT FAMILY                    #
! #                                                         #
! #      (LInearized Discrete Ordinate Radiative Transfer)  #
! #        --           -            -        -        -    #
! #                                                         #
! ###########################################################

! ###########################################################
! #                                                         #
! #  Author :      Robert. J. D. Spurr                      #
! #                                                         #
! #  Address :     RT Solutions, Inc.                       #
! #                9 Channing Street                        #
! #                Cambridge, MA 02138, USA                 #
! #                                                         #
! #  Tel:          (617) 492 1183                           #
! #  Email :        rtsolutions@verizon.net                 #
! #                                                         #
! #  This Version :   3.7 F90                               #
! #  Release Date :   June 2014                             #
! #                                                         #
! #       NEW: THERMAL SUPPLEMENT INCLUDED     (3.2)        #
! #       NEW: OUTGOING SPHERICITY CORRECTION  (3.2)        #
! #       NEW: TOTAL COLUMN JACOBIANS          (3.3)        #
! #       VLIDORT COMPATIBILITY                (3.4)        #
! #       THREADED/OPTIMIZED F90 code          (3.5)        #
! #       EXTERNAL SS / NEW I/O STRUCTURES     (3.6)        #
! #                                                         #
! #       Surface-leaving, BRDF Albedo-scaling     (3.7)    #
! #       Taylor series, BBF Jacobians, ThreadSafe (3.7)    #
! #                                                         #
! ###########################################################

!    #####################################################
!    #                                                   #
!    #   This Version of LIDORT comes with a GNU-style   #
!    #   license. Please read the license carefully.     #
!    #                                                   #
!    #####################################################

! ###############################################################
! #                                                             #
! # Subroutines in this Module                                  #
! #                                                             #
! #            LISPARSE_FUNCTION_PLUS                           #
! #            LIDENSE_FUNCTION_PLUS                            #
! #            HAPKE_FUNCTION_PLUS                              #
! #            RAHMAN_FUNCTION_PLUS                             #
! #            COXMUNK_FUNCTION_PLUS                            #
! #            COXMUNK_FUNCTION_MSR_PLUS                        #
! #                                                             #
! # New BPDF Subroutines in this Module (Version 3.7)           #
! #                                                             #
! #            BPDFVEGN_FUNCTION_PLUS                           #
! #            BPDFSOIL_FUNCTION_PLUS                           #
! #            BPDFNDVI_FUNCTION_PLUS                           #
! #            FRESNEL_SCALAR_PLUS (Private, called by BPDF)    #
! #                                                             #
! # New Cox-Munk Subroutines in this Module (Version 3.7)       #
! #                                                             #
! #            BRDF_Generalized_Glint_plus                      #
! #            BRDF_WhiteCap_Reflectance_plus                   #
! #                                                             #
! ###############################################################

      MODULE brdf_LinSup_kernels_m

!  Rob Extension 12/2/14. BPDF Kernels (replace BREONVEG, BREONSOIL)

      use brdf_sup_aux_m, only : derfc_e, BRDF_Fresnel_Complex

      PRIVATE
      PUBLIC :: LISPARSE_FUNCTION_PLUS, &
                LIDENSE_FUNCTION_PLUS, &
                HAPKE_FUNCTION_PLUS, &
                RAHMAN_FUNCTION_PLUS, &
                COXMUNK_FUNCTION_PLUS, &
                COXMUNK_FUNCTION_MSR_PLUS,    &
                BPDFVEGN_FUNCTION_PLUS,       &
                BPDFSOIL_FUNCTION_PLUS,       &
                BPDFNDVI_FUNCTION_PLUS,       &
                BRDF_Generalized_Glint_plus,  &
                BRDF_WhiteCap_Reflectance_plus

      CONTAINS

      SUBROUTINE LISPARSE_FUNCTION_PLUS              &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,   &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,    &
              LISPARSE_KERNEL, LISPARSE_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, HALF, PIE, TWO

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: LISPARSE_KERNEL
      REAL(fpk), intent(out) :: LISPARSE_DERIVATIVES ( MAXPARS )

!  local variables

      INTEGER    :: J
      REAL(fpk)  :: X_INC, X_REF, SX_INC, SX_REF, ANG_P, TX
      REAL(fpk)  :: T_INC, T_REF, T_INC_SQ, T_REF_SQ
      REAL(fpk)  :: CKSI, DELTA, T, COST, SINT, DSQ, SINTCOST
      REAL(fpk)  :: A, B, H, R, P, Q, DT1, DT2, DT2SQ, QR
      REAL(fpk)  :: DX_H, DX_R, DX_Q, DX_P, DX_QR, DY_Q, TERM2       ! Variables A2, R2 removed 10/12/12, TERM2 added
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise
!    -- Return for special case

      LISPARSE_KERNEL      = ZERO
!mick fix 8/8/2014 - initialize all elements
      LISPARSE_DERIVATIVES = ZERO
      !DO J = 1, NPARS
      !  LISPARSE_DERIVATIVES(J) = ZERO
      !ENDDO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Rob Fix 9/25/14. Hot spot condition  was wrongly inserted here
!      IF ( ( XI .EQ. XJ ) .AND. ( CKPHI.EQ.ONE ) ) RETURN

!  Function
!  ========

!  .. incidence

      TX       = SXJ / XJ
      T_INC    = PARS(2) * TX
      T_INC_SQ = T_INC * T_INC
      ANG_P    = DATAN ( T_INC )
      X_INC    = DCOS(ANG_P)
      SX_INC   = DSIN(ANG_P)

!  .. reflection

      TX       = SXI / XI
      T_REF    = PARS(2) * TX
      T_REF_SQ = T_REF * T_REF
      ANG_P    = DATAN ( T_REF )
      X_REF    = DCOS(ANG_P)
      SX_REF   = DSIN(ANG_P)

!  ksi cosine

      TERM2 = SX_INC * SX_REF * CKPHI
      CKSI  = X_INC  * X_REF + TERM2

!  contributions P and R and derivatives (if flagged)

!    Bug found by Huan Ye (BIRA), Fixed by R. Spurr 12 October 2012
!      For the P term, Division by X_INC omitted in earlier version
!      Linearization of P has changed in the new version

!  Old   P = ( ONE + CKSI ) / X_REF
      P = ( ONE + CKSI ) / X_REF / X_INC
      A = ( ONE / X_INC )
      B = ( ONE / X_REF )
      R = A + B

      DX_R = ZERO
      DX_P = ZERO
      IF ( DO_DERIV_PARS(2) ) THEN
        DX_R = R * ( ONE - ( ONE / A / B ) ) / PARS(2)
! Old       R2   = R * R
! Old       A2   = A * A
! Old       DX_P = ( P * ( ONE + A2 ) - ( R2 / B ) ) / PARS(2) / A2
        DX_P = ( ( SX_INC * SX_INC + SX_REF * SX_REF ) * ( P - ONE ) + &
                  TERM2 * ( X_INC * B + X_REF * A ) ) / PARS(2)
      ENDIF

!  evaluate cos(t)

      DT1   = T_REF_SQ + T_INC_SQ
      DT2   = T_INC * T_REF
      DT2SQ = DT2 * DT2
      DELTA = DSQRT ( DT1 - TWO * DT2 * CKPHI )
      DSQ   = DELTA * DELTA
      H     = DSQRT ( DSQ + SKPHI * SKPHI * DT2SQ )
      COST  = PARS(1) * H / R

!  set Q function and its derivatives if flagged

      DX_Q = ZERO
      DY_Q = ZERO
      IF ( COST .GT. ONE ) THEN
        Q = ONE
      ELSE
        T        = DACOS(COST)
        SINT     = DSQRT ( ONE - COST * COST )
        SINTCOST = SINT * COST
        Q = ONE -  ( ( T - SINTCOST ) / PIE )
        IF ( DO_DERIV_PARS(2) ) THEN
          DX_H = ( TWO * H * H - DSQ ) / H / PARS(2)
          DX_Q = TWO * SINTCOST * ( (DX_H/H) - (DX_R/R) ) / PIE
        ENDIF
        IF ( DO_DERIV_PARS(1) ) THEN
          DY_Q = TWO * SINTCOST / PIE / PARS(1)
        ENDIF
      ENDIF

!  set the kernel
!  --------------

      QR = Q * R
      LISPARSE_KERNEL = HALF * P - QR

!  Set derivatives
!  ---------------

      IF ( DO_DERIV_PARS(1) ) THEN
        LISPARSE_DERIVATIVES(1) = - R * DY_Q
      ENDIF

      IF ( DO_DERIV_PARS(2) ) THEN
        DX_QR = DX_R * Q + DX_Q * R
        LISPARSE_DERIVATIVES(2) = HALF * DX_P - DX_QR
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE LISPARSE_FUNCTION_PLUS

!

      SUBROUTINE LIDENSE_FUNCTION_PLUS               &
             ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,  &
               XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,   &
               LIDENSE_KERNEL, LIDENSE_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, PIE, TWO

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: LIDENSE_KERNEL
      REAL(fpk), intent(out) :: LIDENSE_DERIVATIVES ( MAXPARS )

!  local variables

      INTEGER    :: J
      REAL(fpk)  :: X_INC, X_REF, SX_INC, SX_REF, ANG_P, TX
      REAL(fpk)  :: T_INC, T_REF, T_INC_SQ, T_REF_SQ
      REAL(fpk)  :: CKSI, DELTA, T, COST, SINT, DSQ, SINTCOST
      REAL(fpk)  :: A, B, H, R, P, Q, DT1, DT2, DT2SQ, P_QR
      REAL(fpk)  :: DX_H, DX_R, DX_Q, DX_P, DX_P_QR, DY_Q, TERM2       ! Variables A2, R2 removed 10/12/12, TERM2 added
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise
!    -- Return for special case

      LIDENSE_KERNEL      = ZERO
!mick fix 8/8/2014 - initialize all elements
      LIDENSE_DERIVATIVES = ZERO
      !DO J = 1, NPARS
      !  LIDENSE_DERIVATIVES(J) = ZERO
      !ENDDO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Rob Fix 9/25/14. Hot spot condition  was wrongly inserted here
!      IF ( ( XI .EQ. XJ ) .AND. ( CKPHI.EQ.ONE ) ) RETURN

!  Function
!  ========

!  .. incidence

      TX       = SXJ / XJ
      T_INC    = PARS(2) * TX
      T_INC_SQ = T_INC * T_INC
      ANG_P    = DATAN ( T_INC )
      X_INC    = DCOS(ANG_P)
      SX_INC   = DSIN(ANG_P)

!  .. reflection

      TX       = SXI / XI
      T_REF    = PARS(2) * TX
      T_REF_SQ = T_REF * T_REF
      ANG_P    = DATAN ( T_REF )
      X_REF    = DCOS(ANG_P)
      SX_REF   = DSIN(ANG_P)

!  ksi cosine

      TERM2 = SX_INC * SX_REF * CKPHI
      CKSI  = X_INC  * X_REF + TERM2

!  contributions P and R and derivatives (if flagged)

!    Bug found by Huan Ye (BIRA), Fixed by R. Spurr 12 October 2012
!      For the P term, Division by X_INC omitted in earlier version
!      Linearization of P has changed in the new version

!  Old   P = ( ONE + CKSI ) / X_REF
      P = ( ONE + CKSI ) / X_REF / X_INC
      A = ( ONE / X_INC )
      B = ( ONE / X_REF )
      R = A + B

      DX_R = ZERO
      DX_P = ZERO
      IF ( DO_DERIV_PARS(2) ) THEN
        DX_R = R * ( ONE - ( ONE / A / B ) ) / PARS(2)
! Old       R2   = R * R
! Old       A2   = A * A
! Old       DX_P = ( P * ( ONE + A2 ) - ( R2 / B ) ) / PARS(2) / A2
        DX_P = ( ( SX_INC * SX_INC + SX_REF * SX_REF ) * ( P - ONE ) + &
                  TERM2 * ( X_INC * B + X_REF * A ) ) / PARS(2)
      ENDIF

!  evaluate cos(t)

      DT1   = T_REF_SQ + T_INC_SQ
      DT2   = T_INC * T_REF
      DT2SQ = DT2 * DT2
      DELTA = DSQRT ( DT1 - TWO * DT2 * CKPHI )
      DSQ   = DELTA * DELTA
      H     = DSQRT ( DSQ + SKPHI * SKPHI * DT2SQ )
      COST  = PARS(1) * H / R

!  set Q function and its derivatives if flagged

      DX_Q = ZERO
      DY_Q = ZERO
      IF ( COST .GT. ONE ) THEN
        Q = ONE
      ELSE
        T        = DACOS(COST)
        SINT     = DSQRT ( ONE - COST * COST )
        SINTCOST = SINT * COST
        Q = ONE -  ( ( T - SINTCOST ) / PIE )
        IF ( DO_DERIV_PARS(2) ) THEN
          DX_H = ( TWO * H * H - DSQ ) / H / PARS(2)
          DX_Q = TWO * SINTCOST * ( (DX_H/H) - (DX_R/R) ) / PIE
        ENDIF
        IF ( DO_DERIV_PARS(1) ) THEN
          DY_Q = TWO * SINTCOST / PIE / PARS(1)
        ENDIF
      ENDIF

!  set the kernel
!  --------------

      P_QR = P / Q / R
      LIDENSE_KERNEL = P_QR - TWO

!  Set derivatives
!  ---------------

      IF ( DO_DERIV_PARS(1) ) THEN
        LIDENSE_DERIVATIVES(1) = - P_QR * DY_Q / Q
      ENDIF

      IF ( DO_DERIV_PARS(2) ) THEN
        DX_P_QR = ( DX_P / P ) - ( DX_R / R ) - ( DX_Q / Q )
        LIDENSE_DERIVATIVES(2) = P_QR * DX_P_QR
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE LIDENSE_FUNCTION_PLUS

!

      SUBROUTINE HAPKE_FUNCTION_PLUS                 &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,   &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,    &
              HAPKE_KERNEL, HAPKE_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, TWO, HALF, QUARTER, PIE

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: HAPKE_KERNEL
      REAL(fpk), intent(out) :: HAPKE_DERIVATIVES ( MAXPARS )

!  Hapke Kernel function.
!    - New version, Fresh Coding
!    - Old version uses DISORT code; for validation.

!  input variables:

!    XI, SXI  : Cosine/Sine of angle of reflection (positive)
!    XJ, SXJ  : Cosine/Sine of angle of incidence (positive)
!    XPHI     : Difference of azimuth angles of incidence and reflection
!    PARS(1)  : single scattering albedo in Hapke's BDR model
!    PARS(2)  : angular width parameter of opposition effect in Hapke's model
!    PARS(3)  : Empirical hot spot multiplier

!  local variables
!    B0_EMPIR : empirical factor to account for the finite size of
!               particles in Hapke's BDR model
!    B_HOT    : term that accounts for the opposition effect
!               (retroreflectance, hot spot) in Hapke's BDR model
!    CTHETA   : cosine of phase angle in Hapke's BDR model
!    GAMMA    : albedo factor in Hapke's BDR model
!    PHASE    : scattering phase function in Hapke's BDR model
!    THETA  : phase angle (radians); the angle between incidence and
!             reflection directions in Hapke's BDR model

!  Local variables

      INTEGER    :: J
      REAL(fpk)  :: CTHETA, THETA, PHASE
      REAL(fpk)  :: HOTSPOT, B0_EMPIR, HELP_HOT, B_HOT
      REAL(fpk)  :: SSALBEDO, GAMMA, REFLEC, FUNCTION
      REAL(fpk)  :: HELP_J, GHELP_J, TERM_J
      REAL(fpk)  :: HELP_I, GHELP_I, TERM_I
      REAL(fpk)  :: TI_TJ, DT1, DT2
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise

      HAPKE_KERNEL      = ZERO
!mick fix 8/8/2014 - initialize all elements
      HAPKE_DERIVATIVES = ZERO
      !DO J = 1, NPARS
      !  HAPKE_DERIVATIVES(J) = ZERO
      !ENDDO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  geometrical part

!  This is the code that is in DISORT - not right, I think.
!        CTHETA = XI * XJ + DABS(SXI) *  DABS(SXJ) * CKPHI

      CTHETA = XI * XJ + SXI * SXJ * CKPHI
      THETA  = DACOS( CTHETA )
      PHASE  = ONE + HALF * CTHETA

!  hot spot parameterization

      HOTSPOT  = PARS(2)
      B0_EMPIR = PARS(3)
      HELP_HOT = HOTSPOT + DTAN ( HALF * THETA )
      B_HOT    = B0_EMPIR * HOTSPOT / HELP_HOT

!  Albedo parameterization

      SSALBEDO = PARS(1)
      GAMMA    = DSQRT ( ONE - SSALBEDO )
      HELP_J   = TWO * XJ
      GHELP_J  = ( ONE + HELP_J * GAMMA )
      TERM_J   = ( ONE + HELP_J ) / GHELP_J
      HELP_I   = TWO * XI
      GHELP_I  = ( ONE + HELP_I * GAMMA )
      TERM_I   = ( ONE + HELP_I ) / GHELP_I
      TI_TJ    = TERM_J * TERM_I

!  Function

      REFLEC       = SSALBEDO * QUARTER / ( XI + XJ )
      FUNCTION     = ( ONE + B_HOT ) * PHASE + TI_TJ - ONE
      HAPKE_KERNEL = REFLEC * FUNCTION

!  ssalbedo derivative

      IF ( DO_DERIV_PARS(1) ) THEN
        DT1 = HAPKE_KERNEL / SSALBEDO
        DT2 = ( HELP_J / GHELP_J ) + ( HELP_I / GHELP_I )
        DT2 = DT2 * TI_TJ * HALF / GAMMA
        HAPKE_DERIVATIVES(1) = DT1 + DT2 * REFLEC
      ENDIF

!  Hotspot  derivative

      IF ( DO_DERIV_PARS(2) ) THEN
        DT1 = B_HOT * ( B0_EMPIR - B_HOT ) / B0_EMPIR / HOTSPOT
        HAPKE_DERIVATIVES(2) = DT1 * REFLEC * PHASE
      ENDIF

!  empirical factor derivative

      IF ( DO_DERIV_PARS(3) ) THEN
        DT1 = B_HOT / B0_EMPIR
        HAPKE_DERIVATIVES(3) = DT1 * REFLEC * PHASE
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE HAPKE_FUNCTION_PLUS

!

      SUBROUTINE RAHMAN_FUNCTION_PLUS              &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,  &
              RAHMAN_KERNEL, RAHMAN_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, TWO, ONEP5, PIE

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: RAHMAN_KERNEL
      REAL(fpk), intent(out) :: RAHMAN_DERIVATIVES ( MAXPARS )

!  local variables

      INTEGER    :: J
      REAL(fpk)  :: T_INC, T_REF, DT1, DT2, D_FACT, D_HELPM
      REAL(fpk)  :: CXI, DELTA, K1_SQ, FACT, D_K0, D_K1, D_K2
      REAL(fpk)  :: GEOM, PHASE, RFAC, RFAC1, K0, K1, K2, HELPR
      REAL(fpk)  :: XPHI, CKPHI, HSPOT, UPPER_LIMIT, HELPG, HELPM
      REAL(fpk), PARAMETER  :: SMALL = 1.0d-04

!  Initialise

      RAHMAN_KERNEL      = ZERO
!mick fix 8/8/2014 - initialize all elements
      RAHMAN_DERIVATIVES = ZERO
      !DO J = 1, NPARS
      !  RAHMAN_DERIVATIVES(J) = ZERO
      !ENDDO
      XPHI  = PIE - PHI
      CKPHI = - CPHI
      IF ( XI.EQ.ZERO .OR. XJ.EQ.ZERO ) RETURN

!  parameters

      K0 = PARS(1)
      K1 = PARS(2)
      K2 = PARS(3)

!  Hot Spot
!  --------

!  Value of hot spot

      FACT = K0 * ( 2.0d0 - K0 )
      FACT = FACT * ( 1.0d0 - K1 ) / ( 1.0d0 + K1 ) / ( 1.0d0 + K1 )
      GEOM = ( 2.0d0 * XJ * XJ * XJ ) ** ( K2 - 1.0d0 )
      HSPOT = FACT * GEOM

!  Upper limit ( 5 times hotspot value ). Follwing comments inserted.
!     This function needs more checking; some constraints are
!     required to avoid albedos larger than 1; in particular,
!     the BDREF is limited to 5 times the hotspot value to
!     avoid extremely large values at low polar angles

      UPPER_LIMIT = 5.0d0 * HSPOT

!  hot spot value

      IF ( DABS(PHI) .LT. SMALL .AND. XI.EQ.XJ ) THEN
        RAHMAN_KERNEL = HSPOT
        RETURN
      ENDIF

!  Use upper limit value at edges (low incidence or reflection)

      IF ( XI.LT.SMALL .OR. XJ.LT.SMALL ) THEN
        RAHMAN_KERNEL = UPPER_LIMIT
        RETURN
      ENDIF

!  Main section
!  ------------

!  geometrical angle xi

      CXI = XI * XJ + SXI * SXJ * CKPHI
      IF ( CXI .GT. ONE ) CXI = ONE

!  Phase function

      K1_SQ = K1 * K1
      HELPM = ONE - K1_SQ
      FACT  = ONE + K1_SQ + TWO * K1 * CXI
      PHASE = HELPM / ( FACT ** ONEP5 )

!  Delta and R-factor

      T_INC = SXI / XI
      T_REF = SXJ / XJ
      DT1   = T_INC*T_INC + T_REF*T_REF
      DT2   = T_INC * T_REF
      DELTA = DSQRT ( DT1 - TWO * DT2 * CKPHI )
      HELPR = ONE / ( ONE + DELTA )
      RFAC  = ( ONE - K0 ) * HELPR
      RFAC1 = ONE + RFAC

!  Geom factor and kernel

      HELPG = XI * XJ * ( XI + XJ )
      GEOM  = HELPG ** ( K2 - ONE)
      RAHMAN_KERNEL = K0 * PHASE * RFAC1 * GEOM

!  K0 derivative

      IF ( DO_DERIV_PARS(1) ) THEN
        D_K0   = ( ONE / K0 ) - ( HELPR / RFAC1 )
        RAHMAN_DERIVATIVES(1) = RAHMAN_KERNEL * D_K0
      ENDIF

!  Phase function derivative

      IF ( DO_DERIV_PARS(2) ) THEN
        D_FACT  =   TWO * K1 + TWO * CXI
        D_HELPM = - TWO * K1
        D_K1    = ( D_HELPM / HELPM ) - ONEP5 * ( D_FACT / FACT )
        RAHMAN_DERIVATIVES(2) = RAHMAN_KERNEL * D_K1
      ENDIF

!  K2 derivative

      IF ( DO_DERIV_PARS(3) ) THEN
        D_K2 = DLOG ( HELPG )
        RAHMAN_DERIVATIVES(3) = RAHMAN_KERNEL * D_K2
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE RAHMAN_FUNCTION_PLUS

!

      SUBROUTINE COXMUNK_FUNCTION_PLUS              &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,  &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,   &
              COXMUNK_KERNEL, COXMUNK_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, MINUS_ONE, TWO, FOUR, HALF, QUARTER, PIE, PIO2
      USE BRDF_SUP_AUX_M

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: COXMUNK_KERNEL
      REAL(fpk), intent(inout) :: COXMUNK_DERIVATIVES ( MAXPARS )

!  Critical exponent taken out

      REAL(fpk), PARAMETER   ::  CRITEXP = 88.0D0

!  Local variables

      INTEGER    :: J
      REAL(fpk)  :: Z, Z1, Z2, Z2_SQ_M1, H1, H2, RP, RL, XMP
      REAL(fpk)  :: A, B, TA, ARGUMENT, PROB, FAC1, FAC2
      REAL(fpk)  :: XPHI, CKPHI, T1_R, T2_R, DCOT_R, T1_I, T2_I, DCOT_I
      REAL(fpk)  :: S1, S2, S3, XXI, XXJ, T1, T2, DCOT
      REAL(fpk)  :: SHADOWI, SHADOWR, SHADOW
      REAL(fpk)  :: H1H2, H2Z2, TA_SQ, DFAC2, DH1, DH2, DRP, DRL
      REAL(fpk)  :: D_S1, D_S2, D_T1, D_T2
      REAL(fpk)  :: D_SHADOWI, D_SHADOWR, D_SHADOW

!  Shadow variables

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!               Remark on Use of shadow effect
!               ------------------------------
!  Shadow effect is controlled by the third parameter. That is, if
!  PARS(3) not equal to then shadow effect will be included.
!    --- NPARS should always be 3 for this Kernel.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!  Initialise

      COXMUNK_KERNEL      = ZERO
!mick fix 8/8/2014 - initialize all elements
      COXMUNK_DERIVATIVES = ZERO
      !DO J = 1, NPARS
      !  COXMUNK_DERIVATIVES(J) = ZERO
      !ENDDO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Kernel

!  ..Scatter angles

! old   Z = - XI * XJ + SXI * SXJ * CKPHI
! old   IF ( Z .LT. MINUS_ONE) Z = MINUS_ONE
! old   Z1 = DACOS(-Z)
! old   Z2 = DCOS(Z1*HALF)

      Z = XI * XJ + SXI * SXJ * CKPHI
      IF ( Z .GT. ONE) Z = ONE
      Z1 = DACOS(Z)
      Z2 = DCOS(Z1*HALF)

!  .. Fresnel coefficients

      Z2_SQ_M1 = Z2 * Z2 + MINUS_ONE
      H1 = PARS(2) * Z2
      H2 = DSQRT ( PARS(2) + Z2_SQ_M1 )
      H1H2 = H1 + H2
      RP = ( H1 - H2 ) / H1H2
      H2Z2 = Z2 + H2
      RL = ( Z2 - H2 ) / H2Z2
      XMP = HALF * ( RP*RP + RL*RL )

!  Coxmunk Function

      A = TWO * Z2
      B = ( XI + XJ ) / A
      IF ( B .GT. ONE ) B = ONE
      A  = PIO2 - DASIN(B)
      TA = DTAN(A)
      TA_SQ = TA * TA
      ARGUMENT = TA_SQ  / PARS(1)
      IF ( ARGUMENT .LT. CRITEXP ) THEN
        PROB = DEXP ( - ARGUMENT )
        FAC1 = PROB / PARS(1)
        FAC2 = QUARTER / XI / ( B ** FOUR )
        COXMUNK_KERNEL = XMP * FAC1 * FAC2 / XJ
      ENDIF

!  inverse slope-squared derivative

      IF ( DO_DERIV_PARS(1) ) THEN
        IF ( ARGUMENT .LT. CRITEXP ) THEN
          DFAC2 = ( PARS(1) - TA_SQ ) / PARS(1) / PARS(1)
          COXMUNK_DERIVATIVES(1) = - COXMUNK_KERNEL * DFAC2
        ENDIF
      ENDIF

!  square refractive index derivative
!  --This section of code was formerly at the end of routine
!    -- Now moved here before the shadowing option
!    -- otherwise derivative will not get done
!         Bug found by V. Natraj. 02 February 2007.

      IF ( DO_DERIV_PARS(2) ) THEN
        IF ( ARGUMENT .LT. CRITEXP ) THEN
          DH1 = Z2
          DH2 = HALF / H2
          DRP = ( DH1 * ( ONE - RP ) - DH2 * ( ONE + RP ) ) / H1H2
          DRL =  - DH2 * ( ONE + RL ) / H2Z2
          DFAC2 = ( RP*DRP + RL*DRL ) / XMP
          COXMUNK_DERIVATIVES(2) = COXMUNK_KERNEL * DFAC2
        ENDIF
      ENDIF

!  No Shadow code if not flagged

      IF ( PARS(3) .EQ. ZERO ) RETURN

!  Shadow code
!  -----------

      S1 = DSQRT(PARS(1)/PIE)
      S3 = ONE/(DSQRT(PARS(1)))
      S2 = S3*S3

      XXI  = XI*XI
      DCOT_I = XI/DSQRT(ONE-XXI)
      T1_I   = DEXP(-DCOT_I*DCOT_I*S2)
!      T2_I   = DERFC(DCOT_I*S3)
      T2_I   = DERFC_E(DCOT_I*S3)
      SHADOWI = HALF * ( S1*T1_I/DCOT_I - T2_I )

      XXJ  = XJ*XJ
      DCOT_R = XJ/DSQRT(ONE-XXJ)
      T1_R   = DEXP(-DCOT_R*DCOT_R*S2)
!      T2_R   = DERFC(DCOT_R*S3)
      T2_R   = DERFC_E(DCOT_R*S3)
      SHADOWR = HALF * ( S1*T1_R/DCOT_R - T2_R )

      SHADOW = ONE/(ONE+SHADOWI+SHADOWR)
      COXMUNK_KERNEL = COXMUNK_KERNEL * SHADOW

!  Update Scalar derivatives
!  -------------------------

!  add the shadow derivative to inverse slope-squared derivative

      IF ( DO_DERIV_PARS(1) ) THEN
        D_S1 = HALF / PIE / S1
        D_S2 = - S2 * S2
        D_T1 = - T1_I * DCOT_I * DCOT_I * D_S2
        D_T2 = S2 * S2 * DCOT_I * S1 * T1_I
        D_SHADOWI = HALF * ( D_S1*T1_I/DCOT_I + S1*D_T1/DCOT_I - D_T2 )
        D_T1 = - T1_R * DCOT_R * DCOT_R * D_S2
        D_T2 = S2 * S2 * DCOT_R * S1 * T1_R
        D_SHADOWR = HALF * ( D_S1*T1_R/DCOT_R + S1*D_T1/DCOT_R - D_T2 )
        D_SHADOW = - SHADOW * SHADOW * ( D_SHADOWI + D_SHADOWR )
        COXMUNK_DERIVATIVES(1) = COXMUNK_DERIVATIVES(1) * SHADOW + &
                                 COXMUNK_KERNEL * D_SHADOW / SHADOW
      ENDIF

!  Refractive index derivative, update

      IF ( DO_DERIV_PARS(2) ) THEN
        COXMUNK_DERIVATIVES(2) = COXMUNK_DERIVATIVES(2) * SHADOW
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE COXMUNK_FUNCTION_PLUS

!

      SUBROUTINE COXMUNK_FUNCTION_MSR_PLUS               &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,       &
              ORDER, N_MUQUAD, N_PHIQUAD,                &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,        &
              X_MUQUAD, W_MUQUAD, SX_MUQUAD, WXX_MUQUAD, &
              X_PHIQUAD, W_PHIQUAD,                      &
              COXMUNK_KERNEL, COXMUNK_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, PIE, DEG_TO_RAD, &
                              MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD
      USE BRDF_SUP_AUX_M

      implicit none

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      INTEGER  , intent(in)  :: ORDER
      INTEGER  , intent(in)  :: N_MUQUAD, N_PHIQUAD
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(in)  :: X_MUQUAD   ( MAX_MSRS_MUQUAD )
      REAL(fpk), intent(in)  :: W_MUQUAD   ( MAX_MSRS_MUQUAD )
      REAL(fpk), intent(in)  :: SX_MUQUAD  ( MAX_MSRS_MUQUAD )
      REAL(fpk), intent(in)  :: WXX_MUQUAD ( MAX_MSRS_MUQUAD )
      REAL(fpk), intent(in)  :: X_PHIQUAD  ( MAX_MSRS_PHIQUAD )
      REAL(fpk), intent(in)  :: W_PHIQUAD  ( MAX_MSRS_PHIQUAD )

      REAL(fpk), intent(out) :: COXMUNK_KERNEL
      REAL(fpk), intent(out) :: COXMUNK_DERIVATIVES ( MAXPARS )

!  local variables
!  ---------------

!  help variables

      integer    :: s, n, k, i, i1, N_phiquad_HALF, ni, ki, nr, kr, q
      REAL(fpk)  :: XM, SXM, XMR, SXMR, XMI, SXMI, sum_pr, sumr, sum, w_p
      REAL(fpk)  :: reflec_0, reflec_s, R0Q, D_R0Q, reflec
      REAL(fpk)  :: phi_sub1, cphi_sub1, sphi_sub1
      REAL(fpk)  :: phi_sub2, cphi_sub2, sphi_sub2
      REAL(fpk)  :: d_reflec_0(3), d_reflec_s(3), d_reflec(3)

!  arrays

      REAL(fpk)  :: R0_QUAD_IN   ( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: R0_OUT_QUAD  ( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: D_R0_QUAD_IN ( MAXPARS, MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: D_R0_OUT_QUAD( MAXPARS, MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )

      REAL(fpk)  :: RHOLD   ( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: D_RHOLD ( MAXPARS, MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )

      REAL(fpk)  :: R0_MSRS_QUAD   ( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: D_R0_MSRS_QUAD ( MAXPARS, MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )

!  Safety first zeroing

      REFLEC_0 = ZERO
      COXMUNK_KERNEL = ZERO

      D_REFLEC_0 = ZERO
      COXMUNK_DERIVATIVES = ZERO

!  Single scattering (zero order), Phi is in degrees here!

      CALL COXMUNK_FUNCTION_PLUS                   &
            ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, &
              XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI,  &
              REFLEC_0, D_REFLEC_0 )

!  Higher orders scattering

      REFLEC   = REFLEC_0
      REFLEC_S = ZERO

      D_REFLEC   = D_REFLEC_0
      D_REFLEC_S = ZERO

!  Quadrature output for first order R/T calculations
!        This will be overwritten as the orders increase

      IF ( ORDER.GE.1 ) THEN
        DO K = 1, N_MUQUAD
          XM  = X_MUQUAD(K)
          SXM = SX_MUQUAD(K)
          DO N = 1, N_PHIQUAD
            PHI_SUB1  = X_PHIQUAD(N)
            CPHI_SUB1 = DCOS(PHI_SUB1)
            SPHI_SUB1 = DSIN(PHI_SUB1)
            PHI_SUB2  = PHI*DEG_TO_RAD - X_PHIQUAD(N)
            CPHI_SUB2 = DCOS(PHI_SUB2)
            SPHI_SUB2 = DSIN(PHI_SUB2)
            CALL COXMUNK_FUNCTION_PLUS                              &
               ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,               &
                 XM, SXM, XI, SXI, PHI_SUB2, CPHI_SUB2, SPHI_SUB2,  &
                 R0_OUT_QUAD(K,N), D_R0_OUT_QUAD(:,K,N) )
            CALL COXMUNK_FUNCTION_PLUS                              &
               ( MAXPARS, NPARS, PARS, DO_DERIV_PARS,               &
                 XJ, SXJ, XM, SXM, PHI_SUB1, CPHI_SUB1, SPHI_SUB1,  &
                 R0_QUAD_IN(K,N),  D_R0_QUAD_IN(:,K,N) )
          ENDDO
        ENDDO
      ENDIF

!  Compute the successive orders of scattering.
!  Compute higher orders.

      DO S = 1, ORDER

!  Compute result for this order

         SUMR = ZERO
         DO K = 1, N_MUQUAD
            SUM_PR = ZERO
            DO N = 1, N_PHIQUAD
               W_P = W_PHIQUAD(N)
               SUM = R0_QUAD_IN(K,N) * R0_OUT_QUAD(K,N)
               SUM_PR = SUM_PR + W_P * SUM
            ENDDO
            SUMR = SUMR + SUM_PR * WXX_MUQUAD(K)
         ENDDO
         REFLEC_S = SUMR

!  Derivatives

         DO Q = 1, 2
            IF ( DO_DERIV_PARS(Q) ) THEN
               SUMR = ZERO
               DO K = 1, N_MUQUAD
                  SUM_PR = ZERO
                  DO N = 1, N_PHIQUAD
                     W_P  = W_PHIQUAD(N)
                     SUM_PR = SUM_PR + W_P *   R0_QUAD_IN(K,N)   * D_R0_OUT_QUAD(Q,K,N) &
                                     + W_P * D_R0_QUAD_IN(Q,K,N) *   R0_OUT_QUAD(K,N)
                  ENDDO
                  SUMR = SUMR + SUM_PR * WXX_MUQUAD(K)
               ENDDO
               D_REFLEC_S(Q) = SUMR
            ENDIF
         ENDDO

!  Finish if reached the scattering order desired

         IF ( S.EQ.ORDER ) GO TO 67

!  Compute Reflectance for next order and update
!    Quad-Quad results get computed each time: very wasteful.
!    Have to do this, as the memory is a killer.

         DO KR = 1, N_MUQUAD
            XMR  = X_MUQUAD(KR)
            SXMR = SX_MUQUAD(KR)
            DO NR = 1, N_PHIQUAD

!  Quad-quad calculations

               DO KI = 1, N_MUQUAD
                  XMI  = X_MUQUAD(KI)
                  SXMI = SX_MUQUAD(KI)
                  DO NI = 1, N_PHIQUAD
                     PHI_SUB1  = X_PHIQUAD(NR) - X_PHIQUAD(NI)
                     CPHI_SUB1 = DCOS(PHI_SUB1)
                     SPHI_SUB1 = DSIN(PHI_SUB1)
                     CALL COXMUNK_FUNCTION_PLUS &
                       ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, &
                         XMI, SXMI, XMR, SXMR, PHI_SUB1, CPHI_SUB1, SPHI_SUB1, &
                         R0_MSRS_QUAD(KI,NI), D_R0_MSRS_QUAD(:,KI,NI) )
                  ENDDO
               ENDDO

!  Multiple reflection

               SUMR = ZERO
               DO KI = 1, N_MUQUAD
                  SUM_PR = ZERO
                  DO NI = 1, N_PHIQUAD
                     W_P  = W_PHIQUAD(NI)
                     R0Q = R0_MSRS_QUAD(KI,NI)
                     SUM_PR = SUM_PR + W_P * R0_QUAD_IN(KI,NI) * R0Q
                  ENDDO
!                  SUMR = SUMR + SUM_PR * W_MUQUAD(KI)
                  SUMR = SUMR + SUM_PR * WXX_MUQUAD(KI)
               ENDDO
               RHOLD(KR,NR) = SUMR

!  Derivatives

               DO Q = 1, 2
                  IF ( DO_DERIV_PARS(Q) ) THEN
                     SUMR = ZERO
                     DO KI = 1, N_MUQUAD
                        SUM_PR = ZERO
                        DO NI = 1, N_PHIQUAD
                           W_P = W_PHIQUAD(NI)
                           R0Q   = R0_MSRS_QUAD(KI,NI)
                           D_R0Q = D_R0_MSRS_QUAD(Q,KI,NI)
                           SUM_PR = SUM_PR + W_P *   R0_QUAD_IN(KI,NI)   * D_R0Q &
                                           + W_P * D_R0_QUAD_IN(Q,KI,NI) *   R0Q
                        ENDDO
                        SUMR = SUMR + SUM_PR * WXX_MUQUAD(K)
                     ENDDO
                     D_RHOLD(Q,KR,NR) = SUMR
                  ENDIF
               ENDDO

!  End KR, NR loops

            ENDDO
         ENDDO

!  Update

         DO KR = 1, N_MUQUAD
            DO NR = 1, N_PHIQUAD
               R0_QUAD_IN(KR,NR) = RHOLD(KR,NR)
            ENDDO
         ENDDO
         DO Q = 1, 2
            IF ( DO_DERIV_PARS(Q) ) THEN
               DO KR = 1, N_MUQUAD
                  DO NR = 1, N_PHIQUAD
                     D_R0_QUAD_IN(Q,KR,NR) = D_RHOLD(Q,KR,NR)
                  ENDDO
               ENDDO
            ENDIF
         ENDDO

!  Continuation point for finishing MSR

 67      CONTINUE

!  Add to total

          REFLEC   = REFLEC + REFLEC_S
          D_REFLEC = D_REFLEC + D_REFLEC_S

!  End scattering order loop

      ENDDO

!  Compute total

      COXMUNK_KERNEL = REFLEC
      DO Q = 1, 2
         IF ( DO_DERIV_PARS(Q) ) THEN
            COXMUNK_DERIVATIVES(Q) = D_REFLEC(Q)
         ENDIF
      ENDDO

!  debug

!      write(34,'(1p4e14.5)') coxmunk_kernel, &
!        dacos(xi)/deg_to_rad, dacos(xj)/deg_to_rad, phi

!  Finish

      RETURN
      END SUBROUTINE COXMUNK_FUNCTION_MSR_PLUS

!

      SUBROUTINE BPDFSOIL_FUNCTION_PLUS  &
   ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, &
     BPDFSOIL_KERNEL, BPDFSOIL_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, TWO, HALF, PIE

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: BPDFSOIL_KERNEL
      REAL(fpk), intent(out) :: BPDFSOIL_DERIVATIVES ( MAXPARS )

!  Local variables

      REAL(fpk)  :: REFSQ, Z, Z1, Z2,  FP, XPHI, CKPHI, ATTEN, FP1, DFP
      REAL(fpk)  :: sgamma, cgamma, calpha, calpha_sq, salpha

!  F-.M. Breon BPDF Soil model (2009).
!   This is just the (1,1) component - no polarization

!  Initialise

      BPDFSOIL_KERNEL      = zero
      BPDFSOIL_DERIVATIVES = zero
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Scatter angles, Fresnel reflection
!      PARS(1) = refractive index squared

      Z = XI * XJ + SXI * SXJ * CKPHI
      IF ( Z .GT. ONE) Z = ONE
      Z1 = ACOS(Z)
      Z2 = COS(Z1*HALF)
      REFSQ = PARS(1) * PARS(1)
      CALL FRESNEL_SCALAR_PLUS ( REFSQ, Z2, FP, DFP ) ; DFP = TWO * PARS(1) * DFP

!  Breon and Mick code
!    Scattering angle (=> gamma = scattering angle/2)
!    Note: 0.5 factor applied in alpha & Fp below
!      scat_angle = DACOS(mus*muv + DSQRT((1._fp_kind - mus*mus) &
!                                 *(1._fp_kind - muv*muv)) &
!                                 *DCOS(phi))
!      gamma = scat_angle/2._fp_kind
!      Z2 = dcos(gamma)

!   Angle of the surface that generates specular reflection from
!  sun to view directions (theta)
!      alpha = DACOS(HALF*(mus+muv)/dcos(gamma))

      calpha    = HALF * (xi + xj) / Z2
      calpha_sq = calpha*calpha
      salpha    = sqrt(one - calpha_sq)

      Fp1 = 0.25_fpk / xi / xj

!  BRDF  with attenuation factor

      cgamma = Z2
      sgamma = sqrt ( one - cgamma * cgamma )
      atten  = one - sgamma
      BPDFSOIL_KERNEL = Fp1 * FP * atten

!  Derivatives

      IF ( DO_DERIV_PARS(1) )  BPDFSOIL_DERIVATIVES(1) = Fp1 * DFP * atten

!     Finish

      RETURN
      END SUBROUTINE BPDFSOIL_FUNCTION_PLUS

!

      SUBROUTINE BPDFVEGN_FUNCTION_PLUS  &
   ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, &
     BPDFVEGN_KERNEL, BPDFVEGN_DERIVATIVES )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, TWO, HALF, PIE

      implicit none

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: BPDFVEGN_KERNEL
      REAL(fpk), intent(out) :: BPDFVEGN_DERIVATIVES ( MAXPARS )

!  Local variables

      REAL(fpk)  :: REFSQ, Z, Z1, Z2, FP, FP0, XPHI, CKPHI, ATTEN, PROJECTIONS
      REAL(fpk)  :: sgamma, cgamma, calpha, calpha_sq, salpha
      REAL(fpk)  :: PLEAF, GS, GV, DFP

!  Data coefficients

      REAL(fpk)  :: PLAGIOPHILE_COEFFS(4)
      DATA PLAGIOPHILE_COEFFS /0.43181098_fpk,  0.011187479_fpk, &
                               0.043329567_fpk, 0.19262991_fpk/

!  F-.M. Breon BPDF vegetation model (2009).
!   This is just the (1,1) component - no polarization

!  Initialize

      BPDFVEGN_KERNEL      = zero
      BPDFVEGN_DERIVATIVES = zero
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Scatter angles, Fresnel reflection
!      PARS(1) = refractive index squared

      Z = XI * XJ + SXI * SXJ * CKPHI
      IF ( Z .GT. ONE) Z = ONE
      Z1 = ACOS(Z)
      Z2 = COS(Z1*HALF)
      REFSQ = PARS(1) * PARS(1)
      CALL FRESNEL_SCALAR_PLUS ( REFSQ, Z2, FP, DFP ) ; DFP = TWO * PARS(1) * DFP

!  Breon and Mick code
!    Scattering angle (=> gamma = scattering angle/2)
!    Note: 0.5 factor applied in alpha & Fp below
!      scat_angle = DACOS(mus*muv + DSQRT((1._fp_kind - mus*mus) &
!                                 *(1._fp_kind - muv*muv)) &
!                                 *DCOS(phi))
!      gamma = scat_angle/2._fp_kind
!      Z2 = dcos(gamma)

!   Angle of the surface that generates specular reflection from
!  sun to view directions (theta)

!      alpha = DACOS(HALF*(mus+muv)/dcos(gamma))

      calpha    = HALF * (xi + xj) / Z2
      calpha_sq = calpha*calpha
      salpha    = sqrt(one - calpha_sq)

! Projection of leaf surface to outgoing direction

      gv = PLAGIOPHILE_COEFFS(1) + xi * &
          (PLAGIOPHILE_COEFFS(2) + xi * &
          (PLAGIOPHILE_COEFFS(3) + PLAGIOPHILE_COEFFS(4)*xi))

! Projection of leaf surface to incident direction

      gs = PLAGIOPHILE_COEFFS(1) + xj * &
          (PLAGIOPHILE_COEFFS(2) + xj * &
          (PLAGIOPHILE_COEFFS(3) + PLAGIOPHILE_COEFFS(4)*xj))

! Probability of leaf orientation (plagiophile distr.)

      Pleaf = 16.0_fpk * calpha_sq * salpha  / pie

! Polarization model for vegetation

      PROJECTIONS =  Gv/xi + Gs/xj
      Fp0 = 0.25_fpk * PLEAF / xi / xj / PROJECTIONS

! BRDF  with attenuation factor

      cgamma = Z2
      sgamma = sqrt ( one - cgamma * cgamma )
      atten  = one - sgamma
      BPDFVEGN_KERNEL = Fp0 * FP * atten

!  Derivatives

      IF ( DO_DERIV_PARS(1) )  BPDFVEGN_DERIVATIVES(1) = Fp0 * DFP * atten

!     Finish

      RETURN
      END SUBROUTINE BPDFVEGN_FUNCTION_PLUS

!

      SUBROUTINE BPDFNDVI_FUNCTION_PLUS &
   ( MAXPARS, NPARS, PARS, DO_DERIV_PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, &
     BPDFNDVI_KERNEL, BPDFNDVI_DERIVATIVES )

!  include file of constants

      USE LIDORT_PARS, only : fpk, ZERO, ONE, TWO, HALF, MINUS_ONE, PIE

      IMPLICIT NONE

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: PARS ( MAXPARS )
      LOGICAL  , intent(in)  :: DO_DERIV_PARS ( MAXPARS )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: BPDFNDVI_KERNEL
      REAL(fpk), intent(out) :: BPDFNDVI_DERIVATIVES ( MAXPARS )

!  Local variables

      REAL(fpk)  :: REFSQ, Z, Z1, Z2, FP, XPHI, CKPHI, ATTEN, NDVI, EXPNDVI, C
      REAL(fpk)  :: DFP, DNDVI, DEXPNDVI, KERNEL
      REAL(fpk)  :: sgamma, cgamma

!  F-.M. Breon BPDF NDVI model (2009).
!   This is just the (1,1) component - no polarization

!  Initialise

      BPDFNDVI_KERNEL      = zero
      BPDFNDVI_DERIVATIVES = zero
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Scatter angles, Fresnel reflection
!      PARS(1) = refractive index squared

      Z = XI * XJ + SXI * SXJ * CKPHI
      IF ( Z .GT. ONE) Z = ONE
      Z1 = ACOS(Z)
      Z2 = COS(Z1*HALF)
      REFSQ = PARS(1) * PARS(1)
      CALL FRESNEL_SCALAR_PLUS ( REFSQ, Z2, FP, DFP ) ; DFP = TWO * PARS(1) * DFP

!  PARS(2) = NDVI
!  Exponential of the NDVI ( Out of range values default to zero )

      NDVI = PARS(2) ; DNDVI = ONE
      IF ( NDVI .GT. One .or. NDVI .lt. MINUS_ONE) THEN
        NDVI = zero ; DNDVI = zero
      ENDIF
      EXPNDVI  = EXP ( - NDVI )
      DEXPNDVI = - EXPNDVI * DNDVI

! attenuation factor

      cgamma = Z2
      sgamma = sqrt ( one - cgamma * cgamma )
      atten  = exp ( - sgamma / cgamma )

!  PARS(3) = Scaling Factor

      C = PARS(3)
      KERNEL = 0.25_fpk * atten / ( xi + xj )
      BPDFNDVI_KERNEL = KERNEL * C * FP * EXPNDVI

!  Derivatives

      IF ( DO_DERIV_PARS(1) ) BPDFNDVI_DERIVATIVES(1) = KERNEL * C * DFP * EXPNDVI
      IF ( DO_DERIV_PARS(2) ) BPDFNDVI_DERIVATIVES(2) = KERNEL * C * FP  * DEXPNDVI
      IF ( DO_DERIV_PARS(3) ) BPDFNDVI_DERIVATIVES(3) = KERNEL * FP * EXPNDVI

!  Finish

      RETURN
      END SUBROUTINE BPDFNDVI_FUNCTION_PLUS

!

      SUBROUTINE FRESNEL_SCALAR_PLUS ( PAR1, Z2, FP, DFP )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, MINUS_ONE, ONE, HALF

      implicit none

!  Subroutine arguments

      REAL(fpk), intent(in)  :: PAR1, Z2
      REAL(fpk), intent(out) :: FP, DFP

!  Local variables

      REAL(fpk)  :: Z2_SQ_M1, H1, H2, RP, RL, DH1, DH2, DRP, DRL

!  Code

      Z2_SQ_M1 = Z2 * Z2 + MINUS_ONE
      H1 = PAR1 * Z2
      H2 = SQRT ( PAR1 + Z2_SQ_M1 )
      RP = ( H1 - H2 ) / ( H1 + H2 )
      RL = ( Z2 - H2 ) / ( Z2 + H2 )
      FP = HALF * ( RP*RP + RL*RL )

      DH1 = Z2   ; DH2 = HALF / H2
      DRP = ( ( DH1 - DH2 ) - RP * ( DH1 + DH2 ) ) / ( H1 + H2 )
      DRL =  - DH2 * ( ONE + RL ) / ( Z2 + H2 )
      DFP = RP * DRP + RL * DRL

!  End

      RETURN
      END SUBROUTINE FRESNEL_SCALAR_PLUS

!

SUBROUTINE BRDF_Generalized_Glint_plus &
         ( DO_ISOTROPIC, DO_SHADOW, DO_COEFFS,    &
           REFRAC_R, REFRAC_I, WINDSPEED,         &
           PHI_W, CPHI_W, SPHI_W,                 &
           XJ, SXJ, XI, SXI, PHI, CPHI, SPHI,     &
           SUNGLINT_COEFFS, DSUNGLINT_COEFFS,     &
           SUNGLINT_REFLEC, DSUNGLINT_REFLEC )

      use lidort_pars, only : fpk, zero, one, two, three, four, half, quarter, PIE

      implicit none

!  Subroutine Input arguments
!  --------------------------

!  Flag for using Isotropic Facet distribution

      LOGICAL  , intent(in)    :: DO_ISOTROPIC

!  Flag for including Shadow effect

      LOGICAL  , intent(in)    :: DO_SHADOW

!  Flag for Calculating Cox-Munk Coefficients
!     Only needs to be done once, so intent(inout)

      LOGICAL  , intent(inout) :: DO_COEFFS

!  Real and imaginary parts of refractive index

      real(fpk), intent(in)    :: REFRAC_R
      real(fpk), intent(in)    :: REFRAC_I

!  Windspeed m/s

      real(fpk), intent(in)    :: WINDSPEED

!  Azimuth between Sun and Wind directions. angle in Radians + Cosine/sine

      real(fpk), intent(in)    :: PHI_W, CPHI_W, SPHI_W

!  Incident and reflected ddirections: sines/cosines. Relative azimuth (angle in radians)

      real(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SPHI

!  Subroutine output arguments
!  ---------------------------

!   Glitter reflectance

      real(fpk), intent(out)   :: SUNGLINT_REFLEC
      real(fpk), intent(out)   :: dSUNGLINT_REFLEC

!  Cox-Munk Coefficients. Intent(inout).

      real(fpk), intent(inout) :: SUNGLINT_COEFFS(7)
      real(fpk), intent(inout) :: DSUNGLINT_COEFFS(7)

!  Local arguments
!  ---------------

      real(fpk), PARAMETER :: CRITEXP = 88.0_fpk
      real(fpk), PARAMETER :: six = two * three, twentyfour = six * four

!  Local variables

      real(fpk)  :: B, ZX, ZY, Z, Z1, Z2, XMP
      real(fpk)  :: TILT, TANTILT, TANTILT_SQ, COSTILT
      real(fpk)  :: ARGUMENT, PROB, FAC2, COEFF, VAR, WSigC, WSigU
      real(fpk)  :: XE, XN, XE_sq, XN_sq, XE_sq_1, XN_sq_1
      real(fpk)  :: XPHI, CKPHI, SKPHI, XPHI_W, CKPHI_W, SKPHI_W
      real(fpk)  :: S1, S2, S3, XXI, XXJ, T1, T2, DCOT
      real(fpk)  :: SHADOWI, SHADOWR, SHADOW

      real(fpk)  :: dARGUMENT, dPROB, dCOEFF, dVAR, dWSigC, dWSigU
      real(fpk)  :: AN, AE, dXE, dXN, dXE_sq, dXN_sq, EXPO, dEXPO
      real(fpk)  :: T3, T4, T5, T6, T7, dT3, dT4, dT5, dT6, dT7

!  Initialise output

      SUNGLINT_REFLEC  = ZERO
      dSUNGLINT_REFLEC = ZERO

!  COmpute coefficients, according to 6S formulation

      IF ( DO_COEFFS ) THEN
         SUNGLINT_COEFFS = zero ; DSUNGLINT_COEFFS = zero
         IF ( DO_ISOTROPIC ) THEN
            SUNGLINT_COEFFS(1)  = 0.003_fpk + 0.00512_fpk * WINDSPEED
            DSUNGLINT_COEFFS(1) = 0.00512_fpk
         ELSE
            SUNGLINT_COEFFS(1) = 0.003_fpk + 0.00192_fpk * WINDSPEED ! sigmaC
            SUNGLINT_COEFFS(2) =             0.00316_fpk * WINDSPEED ! sigmaU
            SUNGLINT_COEFFS(3) = 0.010_fpk - 0.00860_fpk * WINDSPEED ! C21
            SUNGLINT_COEFFS(4) = 0.040_fpk - 0.03300_fpk * WINDSPEED ! C03
            SUNGLINT_COEFFS(5) = 0.400_fpk                           ! C40
            SUNGLINT_COEFFS(6) = 0.230_fpk                           ! C04
            SUNGLINT_COEFFS(7) = 0.120_fpk                           ! C22
            DSUNGLINT_COEFFS(1) = 0.00192_fpk
            DSUNGLINT_COEFFS(2) = 0.00316_fpk
            DSUNGLINT_COEFFS(3) = - 0.00860_fpk
            DSUNGLINT_COEFFS(4) = - 0.03300_fpk
         ENDIF
         DO_COEFFS = .false.
      ENDIF

!  Local angles

      XPHI   = PIE - PHI       ! Not used
!     CKPHI  = - CPHI          ! Original, not correct.

      CKPHI  = + CPHI
      SKPHI  = + SPHI

      XPHI_W  = PHI_W
      CKPHI_W = CPHI_W
      SKPHI_W = SPHI_W

!  Tilt angle

      B  = ONE / ( XI + XJ )
      ZX = - SXI * SKPHI * B
      ZY = ( SXJ + SXI * CKPHI ) * B
      TANTILT_SQ = ZX * ZX + ZY * ZY
      TANTILT    = SQRT ( TANTILT_SQ )
      TILT       = ATAN(TANTILT)
      COSTILT    = COS(TILT)

!  Scatter angle

      Z = XI * XJ + SXI * SXJ * CKPHI
      IF ( Z .GT. ONE) Z = ONE
      Z1 = ACOS(Z)
      Z2 = COS(Z1*HALF)

!  Fresnel
!  -------

       CALL BRDF_Fresnel_Complex ( REFRAC_R, REFRAC_I, Z2, XMP )

!  Anisotropic
!  -----------

      IF ( .not. DO_ISOTROPIC ) THEN

!  Variance

         WSigC = ONE / Sqrt(SUNGLINT_COEFFS(1)) ; dWSigC = - half * WSigC * WSigC * WSigC * DSUNGLINT_COEFFS(1)
         WSigU = ONE / Sqrt(SUNGLINT_COEFFS(2)) ; dWSigU = - half * WSigU * WSigU * WSigU * DSUNGLINT_COEFFS(2)
         VAR   = WSigC * WSigU * HALF ; dVAR = half * ( dWSigC * WSigU + WSigC * dWSigU )

!  angles

         AE = (  CKPHI_W * ZX + SKPHI_W * ZY )
         AN = ( -SKPHI_W * ZX + CKPHI_W * ZY )
         XE = AE * WSigC ; XE_sq = XE * XE ; XE_sq_1 = xe_sq - one ; dXE = AE * dWSigC ; dXE_sq = two * dXE * XE
         XN = AN * WSigU ; XN_sq = XN * XN ; XN_sq_1 = xn_sq - one ; dXN = AN * dWSigU ; dXN_sq = two * dXN * XN

!  GC Coefficient

         T3  = XE_sq_1 * XN * half
         dT3 = ( XE_sq_1 * dXN + dXE_sq * XN ) * half
         T4  = ( XN_sq - three ) * XN / six
         dT4 = ( ( XN_sq - three ) * dXN + dXN_sq * XN ) / six
         T5  = ( XE_sq * XE_sq - six * XE_sq + three ) / twentyfour
         dT5 = ( two * dXE_sq * XE_sq - six * dXE_sq ) / twentyfour
         T6  = ( XN_sq * XN_sq - six * XN_sq + three ) / twentyfour
         dT6 = ( two * dXN_sq * XN_sq - six * dXN_sq ) / twentyfour
         T7  = XE_sq_1 * XN_sq_1 / four
         dT7 = ( dXE_sq * XN_sq_1 + XE_sq_1 * dXN_sq ) / four

         Coeff  = ONE - SUNGLINT_COEFFS(3) * T3 &
                      - SUNGLINT_COEFFS(4) * T4 &
                      + SUNGLINT_COEFFS(5) * T5 &
                      + SUNGLINT_COEFFS(6) * T6 &
                      + SUNGLINT_COEFFS(7) * T7
         dCoeff =  - dSUNGLINT_COEFFS(3) * T3 - SUNGLINT_COEFFS(3) * dT3 &
                   - dSUNGLINT_COEFFS(4) * T4 - SUNGLINT_COEFFS(4) * dT4 &
                                              + SUNGLINT_COEFFS(5) * dT5 &
                                              + SUNGLINT_COEFFS(6) * dT6 &
                                              + SUNGLINT_COEFFS(7) * dT7

!  Probability and finish

         ARGUMENT  = (  XE_sq  +  XN_sq ) * HALF
         dARGUMENT = ( dXE_sq  + dXN_sq ) * HALF
         IF ( ARGUMENT .LT. CRITEXP ) THEN
            EXPO = EXP ( - ARGUMENT ) ; dEXPO = - dARGUMENT * EXPO
            PROB = COEFF * EXPO * VAR ; dPROB =  dCOEFF * EXPO * VAR + COEFF * dEXPO * VAR + COEFF * EXPO * dVAR
            FAC2 = QUARTER / XI / XJ / ( COSTILT ** FOUR )
            SUNGLINT_REFLEC  = XMP * PROB  * FAC2
            dSUNGLINT_REFLEC = XMP * dPROB * FAC2
         ENDIF

      ENDIF

!  Isotropic
!  ---------

      IF ( DO_ISOTROPIC ) THEN

!  Compute Probability and finish

         VAR   = SUNGLINT_COEFFS(1) ; dVAR = dSUNGLINT_COEFFS(1)
         ARGUMENT = TANTILT_SQ / VAR
         dARGUMENT = - ARGUMENT * dVAR / VAR
         IF ( ARGUMENT .LT. CRITEXP ) THEN
            EXPO = EXP ( - ARGUMENT ) ; dEXPO = - dARGUMENT * EXPO
            PROB = EXPO / VAR ; dPROB =  ( dEXPO - PROB * dVAR ) / VAR
            FAC2 = QUARTER / XI / XJ / ( COSTILT ** FOUR )
            SUNGLINT_REFLEC  = XMP * PROB  * FAC2
            dSUNGLINT_REFLEC = XMP * dPROB * FAC2
         ENDIF

      ENDIF

!  No Shadow code if not flagged

      IF ( .not. DO_SHADOW  ) RETURN

!  Shadow code

      S1 = SQRT ( VAR / PIE )
      S3 = ONE / ( SQRT(VAR) )
      S2 = S3 * S3

      XXI  = XI*XI
      DCOT = XI / SQRT ( ONE - XXI )
      T1   = EXP ( - DCOT * DCOT * S2 )
      T2   = DERFC_E ( DCOT * S3 )
!      T2   = DCOT * S3 ; CALL HOMEGROWN_ERRFUNC ( T2 )     !  Error function usage in SLEAVE code
      SHADOWI = HALF * ( S1 * T1 / DCOT - T2 )

      XXJ  = XJ*XJ
      DCOT = XJ / SQRT ( ONE - XXJ )
      T1   = EXP ( - DCOT * DCOT * S2 )
      T2   = DERFC_E ( DCOT * S3 )
!      T2   = DCOT * S3 ; CALL HOMEGROWN_ERRFUNC ( T2 )     !  Error function usage in SLEAVE code
      SHADOWR = HALF * ( S1 * T1 / DCOT - T2 )

      SHADOW = ONE / ( ONE + SHADOWI + SHADOWR )
      SUNGLINT_REFLEC = SUNGLINT_REFLEC * SHADOW

!     Finish

      RETURN
END SUBROUTINE BRDF_Generalized_Glint_plus

!

subroutine BRDF_WhiteCap_Reflectance_plus &
    ( WindSpeed, Wavelength, &
      WC_Reflectance, WC_Lambertian, DWC_Reflectance, DWC_Lambertian )

   use lidort_pars, only : fpk, zero

!  Stand-alone routine for computing the WhiteCap Reflectance
!   Based on 6S code, as updated by A. Sayer (2011)

!  Linearization with respect to Wind-speed

!   Made compatible with LIDORT SURFACE LEAVING code
!   renamed for BRDF code (useful if BRDF and SLEAVE are operating together)
!   R. Spurr, 23 April 2014, 28 April 2014

   implicit none

!  Inputs
!    (Wind speed in [m/s], Wavelength in Microns)

   real(fpk), intent(in)  :: WindSpeed
   real(fpk), intent(in)  :: Wavelength

!  output

   real(fpk), intent(out) :: WC_Reflectance
   real(fpk), intent(out) :: WC_Lambertian
   real(fpk), intent(out) :: DWC_Reflectance
   real(fpk), intent(out) :: DWC_Lambertian

!  Data
!  ----

!  Single precision

   real :: Effective_WCRef(39)

! effective reflectance of the whitecaps (Koepke, 1984)
! These are the original values - superseded, A Sayer 05 Jul 2011.
!      data Effective_WCRef/ &
!     0.220,0.220,0.220,0.220,0.220,0.220,0.215,0.210,0.200,0.190,&
!     0.175,0.155,0.130,0.080,0.100,0.105,0.100,0.080,0.045,0.055,&
!     0.065,0.060,0.055,0.040,0.000,0.000,0.000,0.000,0.000,0.000,&
!     0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/

! effective reflectance of the whitecaps (Frouin et al, 1996)
! Assume linear trends between the node points they give
! This is the spectral shape

      data Effective_WCRef/ &
     1.000,1.000,1.000,1.000,0.950,0.900,0.700,0.550,0.500,0.450,&
     0.400,0.350,0.300,0.250,0.200,0.150,0.100,0.050,0.000,0.000,&
     0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,&
     0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/

!  Local variables
!  ---------------

!  Single precision in the original code

   integer :: iwl, iref
   real    :: Wlb, DWlb, WLP, Ref(39), wspd, wl, Ref_i, Rwc, DRwc

!  Initialize

   WC_Reflectance = zero
   WC_Lambertian  = zero
   DWC_Lambertian = zero

!  Single precision inputs in the original

   wspd = real(WindSpeed)
   wl   = real(Wavelength)

!  Scale data for value of 0.22 in the midvisible.

   DO iref = 1,39
      Ref(iref) = 0.22 * Effective_WCRef(iref)
   ENDDO

!  COMPUTE WHITECAPS REFLECTANCE (LAMBERTIAN)

   Wlb    = 0.0 ; DWlb = 0.0
   IF (wspd .le. 9.25) THEN
       Wlb  = 0.01*((3.18e-03)*((wspd-3.7)**3.0))
       DWlb = 0.03*((3.18e-03)*((wspd-3.7)**2.0))
   ELSE IF (wspd .gt. 9.25) THEN
       Wlb  = 0.01*((4.82e-04)*((wspd+1.8)**3.0))
       DWlb = 0.03*((4.82e-04)*((wspd+1.8)**2.0))
   END IF

! Original whitecap calculation - superseded, A. Sayer 05 Jul 2011.
!      W=2.95e-06*(wspd**3.52)

!  Find data point, Linearly interpolate

   iwl   = 1+int((wl-0.2)/0.1)
   wlp   = 0.5+(iwl-1)*0.1
   Ref_i = Ref(iwl+1) + ( wl-wlp)/0.1*(Ref(iwl)-Ref(iwl+1))
   Rwc   = Wlb*Ref_i
   DRwc  = DWlb*Ref_i

!  Final values

   WC_Lambertian   = real(Wlb,fpk)
   DWC_Lambertian  = real(DWlb,fpk)
   WC_Reflectance  = real(Rwc,fpk)
   DWC_Reflectance = real(DRwc,fpk)

!  Finish

   return
end subroutine BRDF_WhiteCap_Reflectance_plus

! end Module

      END MODULE brdf_LinSup_kernels_m

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
      subroutine fitp(am,x,px,rx,w,Ntraj,dpx,ddp,dr,ddr,u,ap,ar,err)

      implicit real*8(a-h,o-z)

      integer*4,intent(IN) :: Ntraj
      integer*4,parameter :: nb=4
      real*8 :: px(Ntraj), w(Ntraj),x(Ntraj),rx(ntraj)
      real*8 :: dpx(Ntraj),dr(ntraj),ddr(ntraj)
      real*8 :: f(nb),s(nb,nb),c(nb,1),cr(nb,1),mat(nb,nb),u(ntraj),
     +          ap(ntraj),ar(ntraj),ddp(ntraj)
      real*8 :: err(2)

      c = 0d0
      s = 0d0
      cr = 0d0

      do i=1,Ntraj
        f=(/x(i),x(i)**2,x(i)**3,1d0/)
        c(1,1) = c(1,1)+px(i)*f(1)*w(i)
        c(2,1) = c(2,1)+px(i)*f(2)*w(i)
        c(3,1) = c(3,1)+f(3)*px(i)*w(i)
        c(4,1) = c(4,1)+f(4)*px(i)*w(i)
      enddo

      do i=1,Ntraj
        f=(/x(i),x(i)**2,x(i)**3,1d0/)
        cr(1,1) = cr(1,1)+rx(i)*f(1)*w(i)
        cr(2,1) = cr(2,1)+rx(i)*f(2)*w(i)
        cr(3,1) = cr(3,1)+f(3)*rx(i)*w(i)
        cr(4,1) = cr(4,1)+f(4)*rx(i)*w(i)
      enddo

      do m=1,nb
        do n=1,nb
          do i=1,Ntraj
            f=(/x(i),x(i)**2,x(i)**3,1d0/)
            s(m,n)=w(i)*f(m)*f(n)+s(m,n)
          end do
        end do
      end do

! calculate matrix c(t)
      mat = s
        call DPOSV('U',nb,1,mat,nb,c,nb,INFO)
        if(INFO/=0) then
        print *, "info=",info
        print *, "polynomial fitting : matrix fails"
        stop
        end if

        call DPOSV('U',nb,1,s,nb,cr,nb,INFO)
        if(INFO/=0) then
        print *, "info=",info
        print *, "polynomial fitting : matrix fails"
        stop
        end if

!     check the error of poly-fit

!      call errp(ntraj,px,x,w,c,err_p)

! the momentum operator r=cf
      do i=1,Ntraj
        ap(i) = c(1,1)*x(i)+c(2,1)*x(i)**2+c(3,1)*x(i)**3+c(4,1)
        ar(i) = cr(1,1)*x(i)+cr(2,1)*x(i)**2+cr(3,1)*x(i)**3+cr(4,1)
        dpx(i) = c(1,1)+2d0*c(2,1)*x(i)+c(3,1)*3d0*x(i)**2
        ddp(i) = 2d0*c(2,1)*x(i)+6d0*x(i)*c(3,1)
        dr(i) = cr(1,1) + 2d0*x(i)*cr(2,1) + cr(3,1)*3d0*x(i)**2
        ddr(i) = 2d0*cr(2,1) + 6d0*x(i)*cr(3,1)
        u(i) = -(rx(i)**2 + dr(i))/(2d0*am)
      enddo

      err = 0d0
      do i=1,ntraj
        err(1) = err(1) + (px(i)-ap(i))**2*w(i)
        err(2) = err(2) + (rx(i)-ar(i))**2*w(i)
      enddo
!c     &             py(i)*(c(1,2)+2d0*c(3,2)*x(i)+c(5,2)*y(i))/am2
!c          dpy(i) = px(i)*(c(2,1)+2d0*c(4,1)*y(i)+c(5,1)*x(i))/am1+
!c     &             py(i)*(c(2,2)+2d0*c(4,2)*y(i)+c(5,2)*x(i))/am2
!! calculate quantum potential
!c          qp(i) = (c(1,1)+2d0*c(2,1)*x(i)+c(3,1)*3d0*x(i)**2)/2d0/am1
!c     &            (c(2,2)+2d0*c(4,2)*y(i)+c(5,2)*x(i))/2d0/am2
!c          qfx(i) = -(2d0*c(2,1)+6d0*x(i)*c(3,1))/2d0/am1
!c     qfy(i) = -(c(5,1)/2d0/am1+c(4,2)/am2)

      return
      end subroutine

!     --------------------------------
!     check the error of poly-fit of p
!     --------------------------------
      subroutine errp(nt,p,x,w,c,error)
      implicit real*8(a-h,o-z)

      integer*4, intent(in)  :: nt
      real*8,    intent(in)  :: w(nt),p(nt),x(nt),c(4,1)
      real*8,    intent(out) :: error

      error = 0d0

      do i=1,nt
        pfit = c(1,1)*x(i)+c(2,1)*x(i)**2+c(3,1)*x(i)**3+c(4,1)
        error = error + (p(i)-pfit)**2*w(i)
      enddo

      return
      end subroutine

ext: f


text:
MODULE enumeration_types
use num_types
implicit none
private
public derivCryst, opList, LabelRotationList, RotPermList, maxLabLength, GroupList, permList

integer, parameter :: maxLabLength = 500  ! maximum length of the character string labeling

!!A list of group operations for each layer in the
!!tree. I.e., the list of stabilizer groups
type :: GroupList
   type(permList), pointer :: layer(:) => null()
endtype GroupList

!!Stores a list of permutations (for the atomic sites in
!!deriv. structure)
type :: permList
   integer, pointer :: perms(:,:) => null()
endtype permList

type RotPermList ! List of permutations that are equivalent
   integer, pointer :: perm(:,:) => null() ! First index is the permutation number
   integer, pointer :: RotIndx(:) ! Which rotations in the list fix the superlattice
   integer nL ! Number of permutations in the list (perm is nL x nAtoms)
   real(dp), pointer :: v(:,:,:) ! Lattice vectors that return the rotated d-vectors
                                 ! back into the first cell 3xnDxnOp
endtype RotPermList

type opList  ! A list of rotation operations
   real(dp), pointer :: rot(:,:,:)
   real(dp), pointer :: shift(:,:)
endtype opList

type LabelRotationList ! A list of indices corresponding to label rotations
   integer, pointer :: lr(:)
endtype LabelRotationList

type cryst ! a structure with a "real-space" definition
   real(dp) :: LV(3,3)
   real(dp), pointer :: bas(:,:) ! 3xn array of atom positions
   integer, pointer :: aTyp(:) ! array of n integers indicating the atom type of each atomic basis
endtype cryst

! Not used yet. Not clear that it is useful.
type derivCryst
   integer :: diag(3)   ! diagonal elements of the SNF
   real(dp):: pLat(3,3) ! Parent lattice
   integer, pointer :: HNF(:,:,:)! Hermite Normal Forms for each SNF
   integer, pointer :: L(:,:,:)  ! Left SNF transformation
   integer, pointer :: labeling(:,:)! list of configurations
endtype derivCryst

type derivStruct
   integer :: diag(3)   ! diagonal elements of the SNF
   real(dp):: pLat(3,3) ! Parent lattice
   real(dp), pointer :: dVec(:,:) ! a 3xn list of the parent basis vectors
   integer :: nD ! Number of dVectors
   integer :: HNF(3,3)! Hermite Normal Forms for each SNF
   integer :: L(3,3)  ! Left SNF transformation
   integer n ! index, number of parent cells in the super
   integer, pointer :: labeling(:)! list of labels
   integer, pointer :: conc(:) ! Number of atoms of each type in labeling
endtype derivStruct
ENDMODULE

max_stars_repo_name: braydenbekker/enumlib

max_issues_repo_name: braydenbekker/enumlib

max_forks_repo_name: braydenbekker/enumlib


text:
! character_2: from Rice bug report (8/5/2010)
program foo
   character(kind=1,len=2) C1
   character(len=2,kind=1) C2
   character(len=42,kind=1) C3
   print *,"frontend works."
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE AGC(BUFIN,BUFOUT,SCRBUF,SPEC)
C     AGC CALCULATES AND APPLIES AGC (AUTOMATIC GAIN CONTROL) TO AN ARRAY OF REAL
C VALUES (A SEISMIC TRACE).  AGC IS A TYPE OF SCALING SUCH THAT THE AVERAGE
C ABSOLUTE VALUE OF THE AMPLITUDES WITHIN A WINDOW WILL BE A CERTAIN LEVEL
C AFTER SCALING.
C     THIS AGC PROGRAM DIFFERS FROM OTHER AGCS IN THE WAY IT HANDLES ZERO (OR
C NEAR ZERO) INPUT AMPLITUDES.  IF THE WINDOW CONTAINS TOO MANY ZEROES, THE
C CALCULATED MULTIPLIER FOR THAT WINDOW IS IGNORED AND THE LAST ONE IS USED.
C E.G. IN THE CASE OF MUTES PRIOR TO P-BREAKS, THE FIRST MULTIPLIER WILL BE
C CALCULATED FROM SAY 75 PERCENT LIVE VALUES, THUS PREVENTING THE P-BREAKS FROM
C BEING OVERDRIVEN.
C     ANY TRACE WITH ALL AMPLITIDES BELOW A CERTAIN LEVEL (DEAD) WILL BE SET TO
C ZERO.
C
C   AGC VARIABLES:
C   ARGUMENTS:
C   BUFIN  - THE INPUT ARRAY (REAL)
C   BUFOUT -  THE OUTPUT ARRAY (REAL) - THE INPUT AND OUTPUT ARRAYS MAY BE THE
C             SAME.
C   SCRBUF  - A SCRATCH ARRAY (REAL) - MUST BE INDXET LONG
C    SPEC   - SPECIAL CASE SWITCH
C           <0, PRINT THE AVERAGE ABSOLUTE VALUE FOR EACH WINDOW
C           =0, DO NOTHING SPECIAL
C COMMON VARIABLE:
C   INDXST  - THE INDX OF THE FIRST SAMPLE IN BUFIN TO PROCESS.  BUFOUT WILL
C             HAVE UNTOUCHED DATA VALUES PRIOR TO INDXST.  THIS PARAMETER CAN BE
C             NON-ZERO IN ORDER TO AVOID SCALING THE WATER LAYER (SAY CPU TIME).
C             THIS SHOULD BE 1 IN MOST CASES.
C   INDXET - THE INDEX OF BUFIN OF THE LAST SAMPLE TO PROCESS.  THE NUMBER OF
C            SAMPLES IN THE ARRAY (WHEN INDXET IS 1).
C   OLEVEL - THE OUTPUT AMPLITUDE LEVEL.  THE AVERAGE ABSOLUTE VALUE.
C   CLIP   - THE CLIP LEVEL.  ANY AMPLITUDE (ABS VALUE) EXCEEDING THE CLIP
C            VALUE WILL BE SET TO CLIP.
C   NPTS   - THE NUMBER OF SAMPLES IN EACH WINDOW( THE WINDOW LENGTH IN SAMPLES)
C            REMEMBER THAT THE WINDOW IS SLIDE (ADVANCED) BY ONE SAMPLE
C   NDEAD  - THE NUMBER OF DEAD VALUES TO ALLOW IN A WINDOW IN ORDER TO GET A
C            NEW WINDOW MULTIPLIER.
C   DEAD   - A LEVEL BELOW WHICH AN AMPLITUDE IS CONSIDERED DEAD.
C   MIDPT  - THE MID POINT OF THE WINDOW TO RECEIVE THE MULTIPLIER FOR THE WINDOW
C            A FORWARD LOOKING WINDOW WOULD HAVE MIDPT=1
c   agcpct - The percent AGC, expressed as percent/100.  agcpct
c            multiplies the multiplier.  An agcpct < 1 has the
c            effect of lessening the AGC (makes it "softer").
C
C
C    PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, JANUARY 1979
c  mod 13 April 1991 by pch - add agcpct
c  mod June 96 by s.s - change the way agcpct works
c
      COMMON /AGCCOM/ INDXST,INDXET,OLEVEL,CLIP,NPTS,NDEAD,DEAD,MIDPT,
     &     agcpct
      DIMENSION BUFIN(1),BUFOUT(1),SCRBUF(1)
      LOGICAL FIRST
      double precision agctot, agcmean
      FIRST=.TRUE.
      NZEROS=0
      SUM=0.
      NDONE=0
C****
C****       PUT THE ABSOLUTE VALUES INTO SCRBUF
C****
      DO 100 I = 1, INDXET
      SCRBUF(I)=ABS(BUFIN(I))
  100 CONTINUE
c  S.S  percent agc change.  First, find the average absolute value of nonzero
c    samples for the entire trace
      if( agcpct .ne. 1.) then
        agctot = 0
        nnonzero = 0
        do 110 i=1,indxet
          if(scrbuf(i) .ne. 0.) then
            agctot = agctot + abs(scrbuf(i))
            nnonzero = nnonzero + 1
          endif
110     continue
        agcmean = agctot / nnonzero
c        write (*,*) 'agcmean ',agcmean,' agcpct',agcpct
      endif
C****
C****            FIND THE FIRST WINDOW WITH A GOOD MULTIPLIER
C****
      INDX=INDXST-1
      DO 190 I = 1, NPTS                                                !  DO THE FIRST WINDOW
         J = INDX + I
         IF( SCRBUF(J) .GT. DEAD ) THEN
             SUM = SUM + SCRBUF(J)
         ELSE
             NZEROS = NZEROS + 1
         ENDIF
  190 CONTINUE
  200 CONTINUE
      NDONE = NDONE + 1
      IF(NZEROS.LT.NDEAD) GO TO 300
      J=INDXST+NDONE-1
      K = J + NPTS
      IF(K.GT.INDXET) GO TO 250
      IF( SCRBUF(J).LE.DEAD) THEN
          NZEROS = NZEROS - 1
          IF( nzeros .LT. 0 ) nzeros = 0
      ELSE
          SUM = SUM - SCRBUF(J)
      ENDIF
      IF( SCRBUF(K) .LE. DEAD ) THEN
          NZEROS = NZEROS + 1
          IF( nzeros .GT. npts ) nzeros = npts
      ELSE
          SUM = SUM + SCRBUF(K)
      ENDIF
      GO TO 200
C****
C****           THE TRACE IS DEAD HERE
C****
  250 CONTINUE
      DO 260 I = 1, INDXET                                              ! MOVE THE INPUT WITHOUT CHANGING THE AMPLITUDES
  260 BUFOUT(I)=BUFIN(I)
      PRINT 270
  270 FORMAT(' ***  WARNING  ***  NOT ENOUGH NON-ZERO VALUES IN ORDER',
     *   ' TO AGC, NO AGC APPLIED.')
      RETURN
C****
C****
  300 CONTINUE
      TEMP=FLOAT(NPTS-NZEROS)                                           ! CHECK FOR NUMBER SIZE
      IF( TEMP .LE. 0. ) THEN
          AVE = CLIP
      ELSE
          AVE = SUM / TEMP
      ENDIF
      IF(AVE.EQ.0.) AVE=.00000001
      IF( SPEC .LT. 0 ) PRINT *,' ave=',ave,' ndone=',ndone
C****
C****         MULTIPLY THE CENTRE POINT OF THE WINDOW AND ADVANCE THE WINDOW
C****         ONE SAMPLE AND REPEAT
C****           USE THE FIRST MULTIPLIER FOR ALL POINTS LESS THAN THE FIRST MID
C****      POINT
C****
c S.S.  changing way pctagc works
c      xmult = olevel / ave * agcpct
c      xmult = olevel / ave / (2.-agcpct)
      IF( agcpct .EQ. 1. ) THEN
          xmult = olevel / ave
      ELSE
          xmult = olevel / (ave + ((1.-agcpct) * agcmean))
      ENDIF
      IF( FIRST ) THEN
          FIRST = .FALSE.
          N = INDXST + NDONE + MIDPT - 3
          DO 390 I = 1, N
             BUFOUT(I) = BUFIN(I) * XMULT
             IF( CLIP .NE. 0. ) THEN
                 IF( ABS(BUFOUT(I)) .GT. CLIP )
     *               BUFOUT(I) = SIGN(CLIP,BUFIN(I))
             ENDIF
  390     CONTINUE
      ENDIF
  400 CONTINUE
      I = INDXST + NDONE + MIDPT - 2
      BUFOUT(I) = BUFIN(I) * XMULT
      IF( CLIP .NE. 0.) THEN
          IF( ABS(BUFOUT(I)) .GT. CLIP )
     *        BUFOUT(I) = SIGN(CLIP,BUFIN(I))
      ENDIF
      IF( INDXST+NPTS+NDONE-1 .GT. INDXET ) GOTO 1000
      J=INDXST+NDONE-1
      K=J+NPTS
      IF( SCRBUF(J) .LE. DEAD ) THEN
          NZEROS = NZEROS - 1
          IF( nzeros .LT. 0 ) nzeros = 0
      ELSE
          SUM = SUM - SCRBUF(J)
      ENDIF
      IF( SCRBUF(K) .LE. DEAD ) THEN
          NZEROS = NZEROS + 1
          IF( nzeros .GT. npts ) nzeros = npts
      ELSE
          SUM = SUM + SCRBUF(K)
      ENDIF
      NDONE = NDONE + 1
      IF( NZEROS .LT. NDEAD ) GOTO 300
      GOTO 400
C****
C****            FINISHED FINDING XMULTIPLIERS, NOW FINISH THE REMAINDER OF THE
C****       LAST WINDOW
C****
 1000 CONTINUE
      J=INDXST+NDONE+MIDPT-1
      DO 1010 I=J,INDXET
      BUFOUT(I)=BUFIN(I)*XMULT
      IF(CLIP.EQ.0.) GO TO 1010
      IF(ABS(BUFOUT(I)).GT.CLIP) BUFOUT(I)=SIGN(CLIP,BUFIN(I))
 1010 CONTINUE
      RETURN
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

! XSUM_KERNELS_ASCII
!
! USAGE
!   NPROC bin/xsum_kernels_ascii INPUT_FILE OUTPUT_DIR
!
!
! COMMAND LINE ARGUMENTS
!   INPUT_FILE             - text file containing list of kernel directories
!   OUTPUT_PATH            - directory to which summed kernels are written
!
!
! DESCRIPTION
!   Sums kernels from directories specified in INPUT_FILE.
!   Writes the resulting sums to OUTPUT_DIR.
!
!   INPUT_FILE is a text file containing a list of absolute or relative paths to
!   kernel directories, one directory per line. Each kernel directory must
!   contain the following ASCII files, which are generated by SPECFEM2D when
!   adjoint simulation and ascii kernel options are specified:
!
!     proc000000_rhop_alpha_beta_kernel.dat
!     proc000000_rho_kappa_mu.dat
!
!   This routine is modeled after xsum_kernels in SPECFEM3D, in so far as
!   similarities between 2D and 3D packages allow.

program sum_kernels_ascii

  use postprocess_par, only: MAX_STRING_LEN, MAX_KERNEL_PATHS, MAX_KERNEL_NAMES, &
    IIN

  implicit none

  integer, parameter :: NARGS = 2

  integer :: ipath, npath, iker, nlines, i,j
  double precision, allocatable, dimension(:) :: kernel1, kernel2, kernel3
  double precision, allocatable, dimension(:) :: kernel_sum1, kernel_sum2, kernel_sum3
  double precision, allocatable, dimension(:,:) :: coord
  character(len=MAX_STRING_LEN) :: input_file, output_dir, kernel_name
  character(len=MAX_STRING_LEN) :: kernel_paths(MAX_KERNEL_PATHS)
  character(len=MAX_STRING_LEN) :: filename, arg(2), line
  integer :: ier
  double precision :: dummy1, dummy2

  write(*,*) 'Running XSUM_KERNELS_ASCII'
  write(*,*)

  if (command_argument_count() /= NARGS) then
    print *, 'USAGE: bin/xsum_kernels_ascii INPUT_FILE OUTPUT_DIR'
    print *, ''
    stop 'Please check command line arguments'
  endif

  ! parse command line arguments
  do i = 1, NARGS
    call get_command_argument(i,arg(i), status=ier)
  enddo

  read(arg(1),'(a)') input_file
  read(arg(2),'(a)') output_dir

  ! parse paths from INPUT_FILE
  npath=0
  open(unit = IIN, file = trim(input_file), status = 'old',iostat = ier)
  if (ier /= 0) then
     print *,'Error opening ',trim(input_file)
     stop 'Please check command line argument: INPUT_FILE'
  endif
  do while (1 == 1)
     read(IIN,'(a)',iostat=ier) line
     if (ier /= 0) exit
     npath = npath+1
     if (npath > MAX_KERNEL_PATHS) stop 'Error number of paths exceeds MAX_KERNEL_PATHS'
     kernel_paths(npath) = line
  enddo
  close(IIN)

  ! print status update
  write(*,*) 'Running XSUM_KERNELS_ASCII'
  write(*,*)

  write(*,*) 'Combining arrays from ',npath,' paths'
  write(*,*)

  ! allocate arrays
  call get_number_gll_points((kernel_paths(1)), nlines)
  allocate( coord(2,nlines) )
  allocate( kernel1(nlines) )
  allocate( kernel2(nlines) )
  allocate( kernel3(nlines) )
  allocate( kernel_sum1(nlines) )
  allocate( kernel_sum2(nlines) )
  allocate( kernel_sum3(nlines) )
  kernel1(:) = 0.0d0
  kernel2(:) = 0.0d0
  kernel3(:) = 0.0d0
  kernel_sum1(:) = 0.0d0
  kernel_sum2(:) = 0.0d0
  kernel_sum3(:) = 0.0d0

  ! sum kernels
  do iker = 1, 2

    if (iker == 1) write(kernel_name,'(a)') 'proc000000_rhop_alpha_beta_kernel.dat'
    if (iker == 2) write(kernel_name,'(a)') 'proc000000_rho_kappa_mu_kernel.dat'

    ! sum kernels
    do ipath = 1, npath
       filename = trim(kernel_paths(ipath)) //'/'// trim(kernel_name)
       open(unit=3,file=filename, status='old', action='read')
       do j = 1,nlines
          read(3,*) dummy1, dummy2, kernel1(j), kernel2(j), kernel3(j)
       enddo
       close(3)

       if (ipath == 1) then
          kernel_sum1 = kernel1
          kernel_sum2 = kernel2
          kernel_sum3 = kernel3
       else
         kernel_sum1(:)=kernel_sum1(:)+kernel1(:)
         kernel_sum2(:)=kernel_sum2(:)+kernel2(:)
         kernel_sum3(:)=kernel_sum3(:)+kernel3(:)
       endif
    enddo

    ! save result
    filename = trim(output_dir) //'/'// trim(kernel_name)
    open(unit=4,file=filename,status='unknown',action='write')
    do j = 1, nlines
          write(4,'(5e11.3)') coord(1,j),coord(2,j),kernel_sum1(j),kernel_sum2(j),kernel_sum3(j)
    enddo
    close(4)

  enddo

  write(*,*) 'Finished writing kernels '
  write(*,*)

end program sum_kernels_ascii

! ------------------------------------------------------------------------------

subroutine get_number_gll_points(kernel_path, nlines)

  use postprocess_par, only: MAX_STRING_LEN, MAX_LINES

  implicit none

  double precision :: dummy1, dummy2, dummy3, dummy4, dummy5

  character(len=MAX_STRING_LEN) :: kernel_path, filename

  integer :: j, ios, nlines

  filename = trim(kernel_path)//'/'//'proc000000_rhop_alpha_beta_kernel.dat'
  open(unit=3,file=filename,status='old',action='read')
  nlines = 0
  do j=1,MAX_LINES
     read(3,*,iostat=ios) dummy1, dummy2, dummy3, dummy4, dummy5
     if (ios /= 0) exit
     nlines=nlines+1
  enddo
  close(3)

end subroutine get_number_gll_points


text:
#define hasShadowMap0 %hasShadowMap0%
#define hasShadowMap1 %hasShadowMap1%
#define hasShadowMap2 %hasShadowMap2%
#define hasShadowMap3 %hasShadowMap3%
#define hasShadowMap4 %hasShadowMap4%
#define hasShadowMap5 %hasShadowMap5%
#define hasShadowMap6 %hasShadowMap6%
#define hasShadowMap7 %hasShadowMap7%

uniform mat4 modelViewInv;

%LIGHTLOOP:START%

#if hasShadowMap%LIGHTID%
	uniform sampler2DShadow shadowMap%LIGHTID%;
#endif

%LIGHTLOOP:END%

varying vec4 ambientGlobal;

uniform sampler2D colorBuffer;
uniform sampler2D normalBuffer;
uniform sampler2D positionBuffer;
uniform sampler2D lumiBuffer;

uniform vec3 camPos;

//uniform mat4 viewMatrix;
//uniform mat4 projMatrix;
//
//uniform vec3 buffersize;
//uniform vec3 camerarange;

float shadowVal(vec4 shadowProj, sampler2DShadow shadowMap, float rad, float shadowRes)
{
	vec2 filterTaps[12];
	filterTaps[0] = vec2(-0.326212,-0.40581);
	filterTaps[1] = vec2(-0.840144,-0.07358);
	filterTaps[2] = vec2(-0.695914,0.457137);
	filterTaps[3] = vec2(-0.203345,0.620716);
	filterTaps[4] = vec2(0.96234,-0.194983);
	filterTaps[5] = vec2(0.473434,-0.480026);
	filterTaps[6] = vec2(0.519456,0.767022);
	filterTaps[7] = vec2(0.185461,-0.893124);
	filterTaps[8] = vec2(0.507431,0.064425);
	filterTaps[9] = vec2(0.89642,0.412458) ;
	filterTaps[10] =vec2(-0.32194,-0.932615);
	filterTaps[11] =vec2(-0.791559,-0.59771);

	float mapScale = rad / shadowRes;

	vec3 shadowST = shadowProj.xyz / shadowProj.w;

	vec4 shadowColor = shadow2D(shadowMap, shadowST);

	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3( filterTaps[0].s*mapScale,  filterTaps[0].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3( filterTaps[1].s*mapScale,  filterTaps[1].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3( filterTaps[2].s*mapScale,  filterTaps[2].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(-filterTaps[3].s*mapScale,  filterTaps[3].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(-filterTaps[4].s*mapScale,  filterTaps[4].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(-filterTaps[5].s*mapScale,  filterTaps[5].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[6].s*mapScale,  filterTaps[6].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[7].s*mapScale,  filterTaps[7].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[8].s*mapScale,  filterTaps[8].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[9].s*mapScale,  filterTaps[9].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[10].s*mapScale,  filterTaps[10].s*mapScale, 0));
	shadowColor += shadow2D(shadowMap, shadowST.xyz + vec3(filterTaps[11].s*mapScale,  filterTaps[11].s*mapScale, 0));
	shadowColor = shadowColor / 13.0;

	return (shadowProj.w > 0.000) ? shadowColor.r : 1.0;
}

void main()
{
	vec2 texCoord = gl_TexCoord[0].xy;

	vec4 color = vec4(0,0,0,1);

	vec4 ecPosDepth = texture2D(positionBuffer, texCoord);

	if (ecPosDepth.w == 0.0) discard;

	vec3 ecPos = ecPosDepth.xyz;

	vec4 baseColor = texture2D(colorBuffer, texCoord);

	vec4 litColor = ambientGlobal;

	vec4 normalSpec = texture2D(normalBuffer, texCoord);

	vec4 lumiDepth = texture2D(lumiBuffer, texCoord);

	vec3 n;

	n = normalSpec.xyz;

	float specVal = normalSpec.w;

	vec3 l;
	float d;
	float atten;
//	float spotDot;
	float spotEffect;
	vec3 v;
	vec3 h;
    float nDotL, nDotHV, power;
	vec4 diffuse, specular;

///////////////////////////////

%LIGHTLOOP:START%
    l = gl_LightSource[%LIGHTID%].position.xyz - ecPos;
    d = length(l);

    atten = 1.0 / (gl_LightSource[%LIGHTID%].constantAttenuation +
                         gl_LightSource[%LIGHTID%].linearAttenuation * d +
                         gl_LightSource[%LIGHTID%].quadraticAttenuation * d * d);

    l = normalize(l);

    spotEffect = -dot(-l,normalize(gl_LightSource[%LIGHTID%].spotDirection));

    spotEffect = (spotEffect > gl_LightSource[%LIGHTID%].spotCosCutoff) ? pow(spotEffect, gl_LightSource[%LIGHTID%].spotExponent):0.0;

    atten *= spotEffect;

    v = normalize(-ecPos);
    //h = normalize(l + v);

    nDotL = max(0.0, dot(n, l));
//    nDotH = max(0.0, dot(n, h));

	nDotHV = max(dot(n,normalize(gl_LightSource[%LIGHTID%].halfVector.xyz)),0.0);

		power = (nDotL == 0.0) ? 0.0 : pow(nDotHV, lumiDepth.z);

		diffuse = gl_FrontLightProduct[%LIGHTID%].diffuse * nDotL * atten;

		specular = specVal * gl_FrontLightProduct[%LIGHTID%].specular * power * atten;

    litColor = baseColor*(diffuse+specular);

#if hasShadowMap%LIGHTID%
	vec3 lightDir = normalize(gl_LightSource[%LIGHTID%].spotDirection.xyz);

	vec4 shadowpos = vec4(ecPosDepth.xyz,1.0) - vec4(normalize(ecPosDepth.xyz-camPos)*0.1,0.0) + vec4(normalize(lightDir)*0.1,0.0);

	vec4 shadowpos0 = gl_TextureMatrix[%LIGHTID%]*shadowpos;

 	litColor = litColor*shadowVal(shadowpos0, shadowMap%LIGHTID%, 2.0, 1024.0);
#endif

	color += litColor;

//litColor = vec4(0,0,0,1);
%LIGHTLOOP:END%
///////////////////////////////
	gl_FragColor = color;
	gl_FragDepth = ecPosDepth.w;
}

max_stars_repo_name: cjcliffe/CubicVR

max_issues_repo_name: mcanthony/CubicVR

max_forks_repo_name: mcanthony/CubicVR


text:
Founded (and originally known as) the Southern Branch of the University of California in 1919, the University of California, Los Angeles, more commonly known as simply UCLA, has occupied a 419acre campus in the Westwood section of Los Angeles since 1929. The campus, bounded by the tony Los Angeles neighborhoods of Brentwood, Westwood, Little Holmby Hills, and Bel Air, is renowned for its professional schools in film, medicine, and law; its undergraduate programs in many disciplines are highly rated as well.

Currently, UCLA has 11 colleges and divisions:
  Anderson School of Management
  Graduate School of Education and Information Studies
  The Henry Samueli School of Engineering and Applied Science
  School of the Arts and Architecture
  School of Dentistry
  School of Law
  School of Medicine
  School of Nursing
  School of Public Health
  School of Public Policy and Social Research
  School of Theater, Film and Television

Like other UC campuses, UCLAs 179 buildings are continually under construction/renovation. Some of its more noted buildings include the original four buildings (Royce, Kinsey, Haynes, and Powell), which often appear in films as a unnamed or fictional Ivy Leaguestyle campus. Other buildings include Pauley Pavilion, Kaufman Hall, and Ackerman Union.

In sports, UCLA is a Division I school, and carries on a rivalry with crosstown USC; its football team practices on campus, but plays their home games at the Rose Bowl in Pasadena. Their mascot is the grizzly bear, with costumed mascots known as Joe and Josie Bruin. This campus pride is reflected in Bruin Walk and the bronze Bruin Bear in Bruin Plaza, at the end of Bruin Walk.

UCLA currently has 24,811 undergraduates and 10,814 graduates (Fall 2005). As of June 2003, there are 331,226 UCLA alumni.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
UC Davis students student

Davis Dollars Intern

Hey Ed... No worries.  The aggressive name got a reaction, which is probably why you chose it (as you say, you made it to protest the front page, which is pretty obvious).  I acknowledged as being aggressive and angry, but the real problem was associating it with a cool project like Davis Dollars, which is what I was lamenting.  Im glad to hear you used the account for Davis Dollars out of honest ignorance (which we all suffer from time to time!) rather than just being nasty.  Thanks for the note, and I look forward to seeing where Davis Dollars goes from here. Users/JabberWokky Evan JabberWokky Edwards, jw.dw@timewarp.org


text:
*DECK ELMHES
      SUBROUTINE ELMHES (NM, N, LOW, IGH, A, INT)
C***BEGIN PROLOGUE  ELMHES
C***PURPOSE  Reduce a real general matrix to upper Hessenberg form
C            using stabilized elementary similarity transformations.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4C1B2
C***TYPE      SINGLE PRECISION (ELMHES-S, COMHES-C)
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine is a translation of the ALGOL procedure ELMHES,
C     NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C     Given a REAL GENERAL matrix, this subroutine
C     reduces a submatrix situated in rows and columns
C     LOW through IGH to upper Hessenberg form by
C     stabilized elementary similarity transformations.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameter, A, as declared in the calling program
C          dimension statement.  NM is an INTEGER variable.
C
C        N is the order of the matrix, A.  N is an INTEGER variable.
C          N must be less than or equal to NM.
C
C        LOW and IGH are two INTEGER variables determined by the
C          balancing subroutine  BALANC.  If  BALANC  has not been
C          used, set LOW=1 and IGH equal to the order of the matrix, N.
C
C        A contains the input matrix.  A is a two-dimensional REAL
C          array, dimensioned A(NM,N).
C
C     On OUTPUT
C
C        A contains the upper Hessenberg matrix.  The multipliers which
C          were used in the reduction are stored in the remaining
C          triangle under the Hessenberg matrix.
C
C        INT contains information on the rows and columns interchanged
C          in the reduction.  Only elements LOW through IGH are used.
C          INT is a one-dimensional INTEGER array, dimensioned INT(IGH).
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  ELMHES
C
      INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
      REAL A(NM,*)
      REAL X,Y
      INTEGER INT(*)
C
C***FIRST EXECUTABLE STATEMENT  ELMHES
      LA = IGH - 1
      KP1 = LOW + 1
      IF (LA .LT. KP1) GO TO 200
C
      DO 180 M = KP1, LA
         MM1 = M - 1
         X = 0.0E0
         I = M
C
         DO 100 J = M, IGH
            IF (ABS(A(J,MM1)) .LE. ABS(X)) GO TO 100
            X = A(J,MM1)
            I = J
  100    CONTINUE
C
         INT(M) = I
         IF (I .EQ. M) GO TO 130
C    .......... INTERCHANGE ROWS AND COLUMNS OF A ..........
         DO 110 J = MM1, N
            Y = A(I,J)
            A(I,J) = A(M,J)
            A(M,J) = Y
  110    CONTINUE
C
         DO 120 J = 1, IGH
            Y = A(J,I)
            A(J,I) = A(J,M)
            A(J,M) = Y
  120    CONTINUE
C    .......... END INTERCHANGE ..........
  130    IF (X .EQ. 0.0E0) GO TO 180
         MP1 = M + 1
C
         DO 160 I = MP1, IGH
            Y = A(I,MM1)
            IF (Y .EQ. 0.0E0) GO TO 160
            Y = Y / X
            A(I,MM1) = Y
C
            DO 140 J = M, N
  140       A(I,J) = A(I,J) - Y * A(M,J)
C
            DO 150 J = 1, IGH
  150       A(J,M) = A(J,M) + Y * A(J,I)
C
  160    CONTINUE
C
  180 CONTINUE
C
  200 RETURN
      END