text:
      SUBROUTINE DT_CKMTX(Ipar,X,Scale2,Pd,F2)
      IMPLICIT NONE
      INTEGER Ipar
      REAL owlam , owlam2 , q02 , q1s , q2 , qq , sb , xx
C**********************************************************************
C
C     PDF based on Regge theory, evolved with .... by ....
C
C     input: IPAR     2212   proton (not installed)
C                       45   Pomeron
C                      100   Deuteron
C
C     output: PD(-6:6) x*f(x)  parton distribution functions
C            (PDFLIB convention: d = PD(1), u = PD(2) )
C
C**********************************************************************

      SAVE
      DOUBLE PRECISION X , Scale2 , Pd(-6:6) , cdn , cup , F2

      INCLUDE 'inc/dtflka'

      DIMENSION qq(7)
C
      q2 = SNGL(Scale2)
      q1s = q2
      xx = SNGL(X)
C  QCD lambda for evolution
      owlam = 0.23
      owlam2 = owlam**2
C  Q0**2 for evolution
      q02 = 2.D0
C
C
C  the conventions are : q(1)=x*u, q(2)=x*d, q(3)=q(4)=x*sbar=x*ubar=...
C                        q(6)=x*charm, q(7)=x*gluon
C
      sb = 0.
      IF ( q2.GT.q02 ) sb = LOG(LOG(q2/owlam2)/LOG(q02/owlam2))
      IF ( Ipar.EQ.2212 ) THEN
         CALL DT_CKMTPR(1,0,xx,sb,qq(1))
         CALL DT_CKMTPR(2,0,xx,sb,qq(2))
         CALL DT_CKMTPR(3,0,xx,sb,qq(3))
         CALL DT_CKMTPR(4,0,xx,sb,qq(4))
         CALL DT_CKMTPR(5,0,xx,sb,qq(5))
         CALL DT_CKMTPR(8,0,xx,sb,qq(6))
         CALL DT_CKMTPR(7,0,xx,sb,qq(7))
C     ELSEIF (IPAR.EQ.45) THEN
C       CALL CKMTPO(1,0,XX,SB,QQ(1))
C       CALL CKMTPO(2,0,XX,SB,QQ(2))
C       CALL CKMTPO(3,0,XX,SB,QQ(3))
C       CALL CKMTPO(4,0,XX,SB,QQ(4))
C       CALL CKMTPO(5,0,XX,SB,QQ(5))
C       CALL CKMTPO(8,0,XX,SB,QQ(6))
C       CALL CKMTPO(7,0,XX,SB,QQ(7))
      ELSE IF ( Ipar.EQ.100 ) THEN
         CALL DT_CKMTDE(1,0,xx,sb,qq(1))
         CALL DT_CKMTDE(2,0,xx,sb,qq(2))
         CALL DT_CKMTDE(3,0,xx,sb,qq(3))
         CALL DT_CKMTDE(4,0,xx,sb,qq(4))
         CALL DT_CKMTDE(5,0,xx,sb,qq(5))
         CALL DT_CKMTDE(8,0,xx,sb,qq(6))
         CALL DT_CKMTDE(7,0,xx,sb,qq(7))
      ELSE

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,A,I4,A)') 'CKMTX:   IPAR =' ,
     &        Ipar , ' not implemented!'
         STOP
      END IF
C
      Pd(-6) = 0.D0
      Pd(-5) = 0.D0
      Pd(-4) = DBLE(qq(6))
      Pd(-3) = DBLE(qq(3))
      Pd(-2) = DBLE(qq(4))
      Pd(-1) = DBLE(qq(5))
      Pd(0) = DBLE(qq(7))
      Pd(1) = DBLE(qq(2))
      Pd(2) = DBLE(qq(1))
      Pd(3) = DBLE(qq(3))
      Pd(4) = DBLE(qq(6))
      Pd(5) = 0.D0
      Pd(6) = 0.D0
      IF ( Ipar.EQ.45 ) THEN
         cdn = (Pd(1)-Pd(-1))/2.D0
         cup = (Pd(2)-Pd(-2))/2.D0
         Pd(-1) = Pd(-1) + cdn
         Pd(-2) = Pd(-2) + cup
         Pd(1) = Pd(-1)
         Pd(2) = Pd(-2)
      END IF
      F2 = 4.0D0/9.0D0*(Pd(2)-Pd(3)+2.0D0*Pd(3))
     &     + 1.0D0/9.0D0*(Pd(1)-Pd(3)+2.0D0*Pd(3))
     &     + 1.0D0/9.0D0*(2.0D0*Pd(3)) + 4.0D0/9.0D0*(2.0D0*Pd(4))
      END SUBROUTINE

ext: f

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
program swap_test
    use :: swap_mod, only : swap
    implicit none
    real :: a = 3.1, b = 4.5
    integer :: i = 13, j = 11
    character(len=20), parameter :: real_fmt = '(A, 2F7.1)'
    character(len=20), parameter :: int_fmt = '(A, 2I3)'

    ! swap on real values
    print real_fmt, 'before', a, b
    call swap(a, b)
    print real_fmt, 'after', a, b

    ! swap on integer values
    print int_fmt, 'before', i, j
    call swap(i, j)
    print int_fmt, 'after', i, j
end program swap_test

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
      PROGRAM FCOORD1
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C INVOKE DEMO DRIVER
C
      CALL COORD1(IERR)
C
C     DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
      STOP
      END
      SUBROUTINE COORD1 (IERROR)
C
C  This example shows the mapping from GKS world coordinates to
C   GKS normalized device coordinates.
C
      DIMENSION XCRA(101),YCRA(101)
      DATA DSL,DSR,DSB,DST / .59,.99,.40,.70 /
      DATA PFL,PFR,PFB,PFT / .71,.98,.42,.69 /
      CALL GSCLIP (0)
      CALL GSLWSC (2.)
C
C Employ the new high quality filled fonts in PLOTCHAR
C
      CALL PCSETC('FN','times-roman')
      CALL PCSETC('FC','%')
C
      CALL SET    (0.,1.,0.,1.,0.,1.,0.,1.,1)
C
      CALL PLCHHQ (.50,.97,
     +     'Mapping from a window in the world coordinate system',
     +     .02,0.,0.)
      CALL PLCHHQ (.50,.93,
     +     'to a viewport in the normalized device coordinate system',
     +     .02,0.,0.)
      CALL PCSETI('FN',21)
      CALL PLCHHQ (.50,.89,
     +'(Used in calls to GKS routines like GCA, GFA, GPL, GPM, and GTX)'
     +,.016,0.,0.)
      CALL PCSETI('FN',29)
      CALL PLCHHQ (.50,.19,
     +     'The window is:  100., 1000., 100., 1000.',
     +     .015,0.,0.)
      CALL PLCHHQ (.50,.15,
     +     'The viewport is:  .15, .95, .10, .90',
     +     .015,0.,0.)
      CALL SET    (.01,.55,.29,.83,0.,1100.,0.,1100.,1)
      CALL PCSETI('FN',22)
      CALL LINE   (   0.,   0.,1100.,   0.)
      CALL LINE   (1100.,   0.,1050.,  25.)
      CALL LINE   (1100.,   0.,1050., -25.)
      CALL PLCHHQ (CFUX(CUFX(1100.)+.015),0.,'X',.015,0.,0.)
      CALL LINE   (   0.,   0.,   0.,1100.)
      CALL LINE   (   0.,1100.,  25.,1050.)
      CALL LINE   (   0.,1100., -25.,1050.)
      CALL PLCHHQ (0.,CFUY(CUFY(1100.)+.015),'Y',.015,0.,0.)
      CALL LINE   ( 100., 100.,1000., 100.)
      CALL LINE   (1000., 100.,1000.,1000.)
      CALL LINE   (1000.,1000., 100.,1000.)
      CALL LINE   ( 100.,1000., 100., 100.)
      CALL PLCHHQ (CFUX(CUFX(100.)+.005),CFUY(CUFY(1000.)-.02),
     +     'WINDOW',.01,0.,-1.)
      CALL PLCHHQ (100.,CFUY(CUFY(100.)-.01),'100',.009,0.,0.)
      CALL PLCHHQ (1000.,CFUY(CUFY(100.)-.01),'1000',.009,0.,0.)
      CALL PLCHHQ (CFUX(CUFX(100.)-.005),100.,'100',.009,0.,1.)
      CALL PLCHHQ (CFUX(CUFX(100.)-.005),1000.,'1000',.009,0.,1.)
      DO 101 I=1,101
         XCRA(I)=200.+7.*REAL(I-1)
         YCRA(I)=550.+425.*SIN(EXP((XCRA(I)-200.)/300.)-1.)
 101  CONTINUE
      CALL CURVE  (XCRA,YCRA,101)
      XCLW=CUFX(100.)
      XCRW=CUFX(1000.)
      YCBW=CUFY(100.)
      YCTW=CUFY(1000.)
      CALL SET    (DSL,DSR,DSB,DST,0.,1.,0.,1.,1)
      CALL LINE   (0.,0.,1.,0.)
      CALL LINE   (1.,0.,1.,1.)
      CALL LINE   (1.,1.,0.,1.)
      CALL LINE   (0.,1.,0.,0.)
      CALL PLCHHQ (CFUX(CUFX(0.)+.005),CFUY(CUFY(1.)-.01),
     +     'DEVICE',.01,0.,-1.)
      CALL SET    (PFL,PFR,PFB,PFT,0.,1.,0.,1.,1)
      CALL LINE   (0.,0.,1.,0.)
      CALL LINE   (1.,0.,1.,1.)
      CALL LINE   (1.,1.,0.,1.)
      CALL LINE   (0.,1.,0.,0.)
      CALL PLCHHQ (CFUX(CUFX(0.)+.005),CFUY(CUFY(1.)-.01),
     +     'PLOTTER FRAME',.01,0.,-1.)
      CALL PLCHHQ (0.,CFUY(CUFY(0.)-.01),'0',.008,0.,0.)
      CALL PLCHHQ (1.,CFUY(CUFY(0.)-.01),'1',.008,0.,0.)
      CALL PLCHHQ (CFUX(CUFX(0.)-.005),0.,'0',.008,0.,1.)
      CALL PLCHHQ (CFUX(CUFX(0.)-.005),1.,'1',.008,0.,1.)
      VPL=PFL+.15*(PFR-PFL)
      VPR=PFL+.95*(PFR-PFL)
      VPB=PFB+.10*(PFT-PFB)
      VPT=PFB+.90*(PFT-PFB)
      CALL SET    (VPL,VPR,VPB,VPT, 100.,1000., 100.,1000.,1)
      CALL LINE   ( 100., 100.,1000., 100.)
      CALL LINE   (1000., 100.,1000.,1000.)
      CALL LINE   (1000.,1000., 100.,1000.)
      CALL LINE   ( 100.,1000., 100., 100.)
      CALL PLCHHQ (CFUX(CUFX(100.)+.005),CFUY(CUFY(1000.)-.01),
     +     'VIEWPORT',.01,0.,-1.)
      CALL PLCHHQ (100.,CFUY(CUFY(100.)-.01),'.15',.008,0.,0.)
      CALL PLCHHQ (1000.,CFUY(CUFY(100.)-.01),'.95',.008,0.,0.)
      CALL PLCHHQ (CFUX(CUFX(100.)-.005),100.,'.10',.008,0.,1.)
      CALL PLCHHQ (CFUX(CUFX(100.)-.005),1000.,'.90',.008,0.,1.)
      CALL CURVE  (XCRA,YCRA,101)
      XCLV=CUFX(100.)
      XCRV=CUFX(1000.)
      YCBV=CUFY(100.)
      YCTV=CUFY(1000.)
      CALL SET    (0.,1.,0.,1.,0.,1.,0.,1.,1)
      CALL PLOTIT (0,0,0)
      CALL DASHDC ('$''',3,1)
      CALL LINED  (XCLW,YCBW,XCLV,YCBV)
      CALL LINED  (XCRW,YCBW,XCRV,YCBV)
      CALL LINED  (XCRW,YCTW,XCRV,YCTV)
      CALL LINED  (XCLW,YCTW,XCLV,YCTV)
      CALL FRAME
      WRITE (6,1001)
      RETURN
C
 1001 FORMAT ('     COORD1 TEST SUCCESSFUL',24X,
     1     'SEE PLOTS TO VERIFY PERFORMANCE')
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! report.f90 --
!     Main source file for the reporting module:
!     write files with various formats (HTML, RTF and LaTex).
!
!     $Id: report.f90,v 1.1 2008-03-25 05:23:36 arjenmarkus Exp $
!
module reporting
    use reporting_rtf
    use reporting_html
    use reporting_latex

    integer, parameter :: report_rtf   = 1
    integer, parameter :: report_html  = 2
    integer, parameter :: report_latex = 3

    type report_data
        integer                 :: type
        type(report_data_rtf)   :: rtf
        type(report_data_html)  :: html
        type(report_data_latex) :: latex
    end type report_data

contains

! report_open --
!     Open the report file
!
! Arguments:
!     report          Data structure identifying the file
!     filename        Name of the file to open
!     title           Title for the document
!     format          Type of document (RTF, HTML or LaTex;
!                     use one of the parameters)
!     papersize       What paper size to use (A4 or Letter)
!
! Note:
!     Capitalisation of the document format is important
!     The papersize is ignored for HTML files
!
subroutine report_open( report, filename, title, format, papersize )

    type(report_data) :: report
    character(len=*)  :: filename
    character(len=*)  :: title
    integer           :: format
    character(len=*)  :: papersize

    select case( format )
        case( report_rtf )
            report%type = format
            call report_open_rtf( report%rtf, filename, title, papersize )
        case( report_html )
            report%type = format
            call report_open_html( report%html, filename, title, papersize )
        case( report_latex )
            report%type = format
            call report_open_latex( report%latex, filename, title, papersize )
        case default
            report%type = 0
            write(*,*) 'Unknown document format: ', format
    end select
end subroutine report_open

include 'report_generated.f90'

end module reporting


!*==sgebak.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b SGEBAK
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download SGEBAK + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE SGEBAK( JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV,
!                          INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          JOB, SIDE
!       INTEGER            IHI, ILO, INFO, LDV, M, N
!       ..
!       .. Array Arguments ..
!       REAL               V( LDV, * ), SCALE( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> SGEBAK forms the right or left eigenvectors of a real general matrix
!> by backward transformation on the computed eigenvectors of the
!> balanced matrix output by SGEBAL.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] JOB
!> \verbatim
!>          JOB is CHARACTER*1
!>          Specifies the type of backward transformation required:
!>          = 'N': do nothing, return immediately;
!>          = 'P': do backward transformation for permutation only;
!>          = 'S': do backward transformation for scaling only;
!>          = 'B': do backward transformations for both permutation and
!>                 scaling.
!>          JOB must be the same as the argument JOB supplied to SGEBAL.
!> \endverbatim
!>
!> \param[in] SIDE
!> \verbatim
!>          SIDE is CHARACTER*1
!>          = 'R':  V contains right eigenvectors;
!>          = 'L':  V contains left eigenvectors.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The number of rows of the matrix V.  N >= 0.
!> \endverbatim
!>
!> \param[in] ILO
!> \verbatim
!>          ILO is INTEGER
!> \endverbatim
!>
!> \param[in] IHI
!> \verbatim
!>          IHI is INTEGER
!>          The integers ILO and IHI determined by SGEBAL.
!>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
!> \endverbatim
!>
!> \param[in] SCALE
!> \verbatim
!>          SCALE is REAL array, dimension (N)
!>          Details of the permutation and scaling factors, as returned
!>          by SGEBAL.
!> \endverbatim
!>
!> \param[in] M
!> \verbatim
!>          M is INTEGER
!>          The number of columns of the matrix V.  M >= 0.
!> \endverbatim
!>
!> \param[in,out] V
!> \verbatim
!>          V is REAL array, dimension (LDV,M)
!>          On entry, the matrix of right or left eigenvectors to be
!>          transformed, as returned by SHSEIN or STREVC.
!>          On exit, V is overwritten by the transformed eigenvectors.
!> \endverbatim
!>
!> \param[in] LDV
!> \verbatim
!>          LDV is INTEGER
!>          The leading dimension of the array V. LDV >= max(1,N).
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup realGEcomputational
!
!  =====================================================================
      SUBROUTINE SGEBAK(Job,Side,N,Ilo,Ihi,Scale,M,V,Ldv,Info)
      IMPLICIT NONE
!*--SGEBAK133
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Job , Side
      INTEGER Ihi , Ilo , Info , Ldv , M , N
!     ..
!     .. Array Arguments ..
      REAL V(Ldv,*) , Scale(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ONE
      PARAMETER (ONE=1.0E+0)
!     ..
!     .. Local Scalars ..
      LOGICAL leftv , rightv
      INTEGER i , ii , k
      REAL s
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
!     ..
!     .. External Subroutines ..
      EXTERNAL SSCAL , SSWAP , XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MAX , MIN
!     ..
!     .. Executable Statements ..
!
!     Decode and Test the input parameters
!
      rightv = LSAME(Side,'R')
      leftv = LSAME(Side,'L')
!
      Info = 0
      IF ( .NOT.LSAME(Job,'N') .AND. .NOT.LSAME(Job,'P') .AND.          &
     &     .NOT.LSAME(Job,'S') .AND. .NOT.LSAME(Job,'B') ) THEN
         Info = -1
      ELSEIF ( .NOT.rightv .AND. .NOT.leftv ) THEN
         Info = -2
      ELSEIF ( N<0 ) THEN
         Info = -3
      ELSEIF ( Ilo<1 .OR. Ilo>MAX(1,N) ) THEN
         Info = -4
      ELSEIF ( Ihi<MIN(Ilo,N) .OR. Ihi>N ) THEN
         Info = -5
      ELSEIF ( M<0 ) THEN
         Info = -7
      ELSEIF ( Ldv<MAX(1,N) ) THEN
         Info = -9
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('SGEBAK',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 ) RETURN
      IF ( M==0 ) RETURN
      IF ( LSAME(Job,'N') ) RETURN
!
      IF ( Ilo/=Ihi ) THEN
!
!     Backward balance
!
         IF ( LSAME(Job,'S') .OR. LSAME(Job,'B') ) THEN
!
            IF ( rightv ) THEN
               DO i = Ilo , Ihi
                  s = Scale(i)
                  CALL SSCAL(M,s,V(i,1),Ldv)
               ENDDO
            ENDIF
!
            IF ( leftv ) THEN
               DO i = Ilo , Ihi
                  s = ONE/Scale(i)
                  CALL SSCAL(M,s,V(i,1),Ldv)
               ENDDO
            ENDIF
!
         ENDIF
      ENDIF
!
!     Backward permutation
!
!     For  I = ILO-1 step -1 until 1,
!              IHI+1 step 1 until N do --
!
      IF ( LSAME(Job,'P') .OR. LSAME(Job,'B') ) THEN
         IF ( rightv ) THEN
            DO ii = 1 , N
               i = ii
               IF ( i<Ilo .OR. i>Ihi ) THEN
                  IF ( i<Ilo ) i = Ilo - ii
                  k = Scale(i)
                  IF ( k/=i ) CALL SSWAP(M,V(i,1),Ldv,V(k,1),Ldv)
               ENDIF
            ENDDO
         ENDIF
!
         IF ( leftv ) THEN
            DO ii = 1 , N
               i = ii
               IF ( i<Ilo .OR. i>Ihi ) THEN
                  IF ( i<Ilo ) i = Ilo - ii
                  k = Scale(i)
                  IF ( k/=i ) CALL SSWAP(M,V(i,1),Ldv,V(k,1),Ldv)
               ENDIF
            ENDDO
         ENDIF
      ENDIF
!
!
!     End of SGEBAK
!
      END SUBROUTINE SGEBAK


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Log: paint.f,v $
C Revision 1.3  2009/03/25 12:36:46  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.2  2009/01/22 21:34:22  gdsjaar
C There were several inline dbnums common blocks. Replaced with the
C include so they all have the same size with the added variable types.
C
C Added minor support for nodeset and sideset variables.
C
C It can print the count and the names, but that is all
C at this time.
C
C Revision 1.1  1994/04/07 20:06:40  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.3  1993/08/04  15:58:41  gdsjaar
c Performed some optimization to speed up blot - Inlined the locrea function
c
c Revision 1.2  1990/12/14  08:54:31  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE PAINT (VARNP, LENF, NLNKF, LINKF, NXFAC, IXFAC,
     &   XN, YN, ZN, XF, YF, ZF, ISVOK, FMIN, FMAX, *)
C=======================================================================

C   --*** PAINT *** (DETOUR) Paint contours
C   --   Modified by John H. Glick - 10/26/88
C   --   Written by Amy Gilkey - revised 03/14/88
C   --
C   --PAINT paints contour sections in a color sequence.  The contour
C   --algorithm assumes that the elements do not contain internal nodes.
C   --The element interpolation field is approximated by logically drawing
C   --lines from each node to the element center and, thusly, dividing the
C   --element into triangles.  Contour sections are then drawn by connecting
C   --the intersection points of the sub-element edges and the contour
C   --plane.
C   --
C   --The element block status indicates which element blocks are active
C   --for contours.  No contours are drawn in inactive elements.
C   --
C   --Parameters:
C   --   VARNP - IN - the contour function values
C   --   LENF - IN - the cumulative face counts by element block
C   --   NLNKF - IN - the number of nodes per face
C   --   LINKF - IN - the connectivity for all faces
C   --   NXFAC - IN - the number of ordered faces
C   --   IXFAC - IN - the indices of the ordered faces
C   --   XN, YN, ZN - IN - the nodal coordinates
C   --   XF, YF, ZF - IN - the face center coordinates
C   --   ISVOK - IN - ISVOK(i) is true iff the contour variable is defined
C   --      for element block i (always true if nodal variable)
C   --   FMIN, FMAX - IN - the minimum and maximum function value
C   --   * - return statement if the cancel function is active
C   --
C   --Common Variables:
C   --   Uses NUMEL, NELBLK of /DBNUMS/
C   --   Uses IS3DIM of /D3NUMS/
C   --   Uses NCNTR, NOCMIN, NOCMAX of /CNTR/

      include 'dbnums.blk'
      COMMON /D3NUMS/ IS3DIM, NNPSUR, NUMNPF, LLNSET
      LOGICAL IS3DIM
      COMMON /CNTR/   CINTOK, LINCON, NCNTR, CMIN, CMAX, DELC,
     &   CINTV(256), NOCMIN, NOCMAX, LABINC, MAXMIN, MAXMAX
      LOGICAL CINTOK, LINCON, NOCMIN, NOCMAX

      REAL VARNP(*)
      INTEGER LENF(0:NELBLK)
      INTEGER NLNKF(NELBLK)
      INTEGER LINKF(*)
      INTEGER IXFAC(*)
      REAL XN(*), YN(*), ZN(*)
      REAL XF(*), YF(*), ZF(*)
      LOGICAL ISVOK(NELBLK)

      LOGICAL GRABRT

C   --Fill the CINTV array with the contour intervals

      IF (.NOT. CINTOK) THEN
         DO 100 NC = 1, NCNTR+1
            CINTV(NC) = CNTRI (NC)
  100    CONTINUE
      END IF

C   --Set the contour minimum and maximum (in case NOCMIN, etc.)
      CNTRMN = FMIN - ABS (DELC)
      CNTRMX = FMAX + ABS (DELC)

      DO 120 IX = 1, NXFAC
         IFAC = IXFAC(IX)
         IELB = 0
         IXL = IDBLNK (IELB, IFAC, LENF, NLNKF)
         IF (NLNKF(IELB) .LE. 2) GOTO 120

         IF ((.NOT. IS3DIM)  .AND.  (NLNKF(IELB) .EQ. 9)) THEN
            NNPF = 8
         ELSE
            NNPF = NLNKF(IELB)
         ENDIF

         IF (ISVOK(IELB)) THEN

C         --Compute the minimum and maximum values for the face
            CALL CFMAX (VARNP, NLNKF(IELB), LINKF(IXL), FEMIN, FEMAX)

C         --Find limiting contour intervals
              difx = abs(cintv(1) - femax)
              difn = abs(cintv(1) - femin)
              imaxc = 1
              iminc = 1
              do 10 i=2, ncntr+1
                difxx = abs(cintv(i) - femax)
                if (difx .gt. difxx) then
                  difx = difxx
                  imaxc = i
                end if
                difnn = abs(cintv(i) - femin)
                if (difn .gt. difnn) then
                  difn = difnn
                  iminc = i
                end if
 10           continue
            IF (DELC .GE. 0.0) THEN
c$$$               MINC = LOCREA (FEMIN, NCNTR+1, CINTV)
c$$$               MAXC = LOCREA (FEMAX, NCNTR+1, CINTV)
              minc = iminc
              maxc = imaxc
              IF (FEMIN .LT. CINTV(MINC)) MINC = MINC - 1
              IF (FEMAX .LT. CINTV(MAXC)) MAXC = MAXC - 1
            ELSE
c$$$               MINC = LOCREA (FEMAX, NCNTR+1, CINTV)
c$$$               MAXC = LOCREA (FEMIN, NCNTR+1, CINTV)
              minc = imaxc
              maxc = iminc
              IF (FEMAX .GE. CINTV(MINC)) MINC = MINC - 1
              IF (FEMIN .GE. CINTV(MAXC)) MAXC = MAXC - 1
            END IF

            IF (NOCMIN .AND. (MINC .LT. 1)) THEN
               MINC = 1
               IF (MINC .GT. MAXC) MAXC = MINC
            END IF
            IF (NOCMAX .AND. (MAXC .GT. NCNTR)) THEN
               MAXC = NCNTR
               IF (MINC .GT. MAXC) MINC = MAXC
            END IF
            IF ((MAXC .LT. 1) .OR. (MINC .GT. NCNTR)) THEN
               MINC = -1
               MAXC = -1
            END IF

         ELSE
C         --Not selected, paint face black
            MINC = -1
            MAXC = -1
         END IF

C      --Skip this contour if the values are outside the contour range

         IF (MINC .EQ. MAXC) THEN

C         --Face is entirely inside one contour area

            CALL GRCOLR (MINC)
            CALL SOLIDF (NNPF, LINKF(IXL), XN, YN, ZN)

         ELSE

C         --Face has several contour areas

            IF ((MINC .LT. 1) .OR. (MAXC .GT. NCNTR)) THEN

C            --Some part of face is outside contour limits, paint black

               CALL GRCOLR (-1)
               MINC = MAX (MINC, 1)
               MAXC = MIN (MAXC, NCNTR)
            END IF

            DO 110 NC = MINC, MAXC
               CALL GRCOLR (NC)
               IF (DELC .GE. 0.0) THEN
                  CNTR0 = CNTRI (NC)
                  IF (NOCMIN .AND. (NC .EQ. 1)) CNTR0 = CNTRMN
                  CNTR1 = CNTRI (NC+1)
                  IF (NOCMAX .AND. (NC .EQ. NCNTR)) CNTR1 = CNTRMX
               ELSE
                  CNTR1 = CNTRI (NC)
                  IF (NOCMIN .AND. (NC .EQ. 1)) CNTR1 = CNTRMX
                  CNTR0 = CNTRI (NC+1)
                  IF (NOCMAX .AND. (NC .EQ. NCNTR)) CNTR0 = CNTRMN
               END IF

               IF ((FEMIN .GE. CNTR1) .OR. (FEMAX .LT. CNTR0))
     &            print *, femin, cntr1, femax, cntr0

               IF (GRABRT ()) RETURN 1
               CALL PAINTF (CNTR0, CNTR1,
     &            VARNP, NLNKF(IELB), LINKF(IXL),
     &            XN, YN, ZN, XF(IFAC), YF(IFAC), ZF(IFAC))
  110       CONTINUE
         END IF
  120 CONTINUE

      CALL PLTFLU

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
c subroutine outm.f writes out data for goldstein last change  6/6/95
c expanded to write out atmos and sea ice data (Bob 10/5/02)
c
c AY (13/07/04) : modified to output surface fluxes
c                 note the organisation of output file (one field after
c                 the other, no clever interlacing)
c
c AY (23/07/04) : modified to reflect new wind speed/stress field
c                 organisation (i.e. flat 2D fields for u and v points)
c
c AY (12/07/05) : removed surplus argument to function
c
      subroutine outm_surf_ocn_sic(unit,
     :        otemp, osaln,
     :        atemp, ashum, apres,
     :        sich, sica, tice,
     :        windspdxu_atm, windspdyu_atm,
     :        net_sw, net_lw,
     :        albedo_ocn,albedo_sic,
     :        stressxu_ocn, stressyu_ocn, usurf,
     :        fxlho, fxsho, fxswo, fxlwo, evap_net,
     :        fxlha, fxsha, evap_atm,
     :        dthsic, dtareasic)

#include "ocean.cmn"

      integer i, j, unit

      real otemp(imax,jmax), osaln(imax,jmax),
     :     atemp(imax,jmax), ashum(imax,jmax), apres(imax,jmax),
     :     sich(imax,jmax), sica(imax,jmax), tice(imax,jmax),
     :     windspdxu_atm(imax,jmax), windspdyu_atm(imax,jmax),
     :     net_sw(imax,jmax), net_lw(imax,jmax),
     :     albedo_ocn(imax,jmax), albedo_sic(imax,jmax),
     :     stressxu_ocn(imax,jmax), stressyu_ocn(imax,jmax),
     :     fxlho(imax,jmax), fxsho(imax,jmax),
     :     fxswo(imax,jmax), fxlwo(imax,jmax), evap_net(imax,jmax),
     :     fxlha(imax,jmax), fxsha(imax,jmax), evap_atm(imax,jmax),
     :     dthsic(imax,jmax), dtareasic(imax,jmax)
c
      real usurf(imax, jmax)
c
c Ocean : surface temperature
      do j=1,jmax
         do i=1,imax
            write(unit,* )otemp(i,j)
         enddo
      enddo
c
c Ocean : surface salinity
      do j=1,jmax
         do i=1,imax
c
c AY (01/08/04) : correct salinity to PSU
            if (k1(i,j).le.kmax) then
               write(unit,* )osaln(i,j)+saln0
            else
               write(unit,* )osaln(i,j)
            endif
         enddo
      enddo
c
c Atmosphere : lowest level temperature
      do j=1,jmax
         do i=1,imax
            write(unit,* )atemp(i,j)
         enddo
      enddo
c
c Atmosphere : lowest level specific humidity
      do j=1,jmax
         do i=1,imax
            write(unit,* )ashum(i,j)
         enddo
      enddo
c
c Atmosphere : lowest level pressure
      do j=1,jmax
         do i=1,imax
            write(unit,* )apres(i,j)
         enddo
      enddo
c
c Sea-ice : height
      do j=1,jmax
         do i=1,imax
            write(unit,* )sich(i,j)
         enddo
      enddo
c
c Sea-ice : fractional area
      do j=1,jmax
         do i=1,imax
            write(unit,* )sica(i,j)
         enddo
      enddo
c
c Sea-ice : surface temperature
      do j=1,jmax
         do i=1,imax
            write(unit,* )tice(i,j)
         enddo
      enddo
c
c Atmosphere : x wind speed at u point
      do j=1,jmax
         do i=1,imax
            write(unit,* )windspdxu_atm(i,j)
         enddo
      enddo
c
c Atmosphere : y wind speed at u point
      do j=1,jmax
         do i=1,imax
            write(unit,* )windspdyu_atm(i,j)
         enddo
      enddo
c
c Atmosphere : net short-wave flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )net_sw(i,j)
         enddo
      enddo
c
c Atmosphere : net long-wave flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )net_lw(i,j)
         enddo
      enddo
c
c Ocean : albedo
      do j=1,jmax
         do i=1,imax
            write(unit,* )albedo_ocn(i,j)
         enddo
      enddo
c
c Sea-ice : albedo
      do j=1,jmax
         do i=1,imax
            write(unit,* )albedo_sic(i,j)
         enddo
      enddo
c
c AY (23/09/04) : no longer calculated in surf_ocn_sic.F
c Ocean and Sea-ice : average albedo
c     do j=1,jmax
c        do i=1,imax
c           write(unit,* )albedo_net(i,j)
c        enddo
c     enddo
c
c Ocean and Sea-ice : average surface temperature
c     do j=1,jmax
c        do i=1,imax
c           write(unit,* )temp_net(i,j)
c        enddo
c     enddo
c
c Ocean : x wind stress at u point
      do j=1,jmax
         do i=1,imax
            write(unit,* )stressxu_ocn(i,j)
         enddo
      enddo
c
c Ocean : y wind stress at u point
      do j=1,jmax
         do i=1,imax
            write(unit,* )stressyu_ocn(i,j)
         enddo
      enddo
c
c Ocean : surface wind speed at tracer point
      do j=1,jmax
         do i=1,imax
            write(unit,* )usurf(i,j)
         enddo
      enddo
c
c Ocean : latent heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxlho(i,j)
         enddo
      enddo
c
c Ocean : sensible heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxsho(i,j)
         enddo
      enddo
c
c Ocean : short-wave heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxswo(i,j)
         enddo
      enddo
c
c Ocean : long-wave heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxlwo(i,j)
         enddo
      enddo
c
c Ocean : evaporation flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )evap_net(i,j)
         enddo
      enddo
c
c Atmosphere : latent heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxlha(i,j)
         enddo
      enddo
c
c Atmosphere : sensible heat flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )fxsha(i,j)
         enddo
      enddo
c
c Atmosphere : evaporation flux
      do j=1,jmax
         do i=1,imax
            write(unit,* )evap_atm(i,j)
         enddo
      enddo
c
c Sea-ice : change in sea-ice height
      do j=1,jmax
         do i=1,imax
            write(unit,* )dthsic(i,j)
         enddo
      enddo
c
c Sea-ice : change in sea-ice area
      do j=1,jmax
         do i=1,imax
            write(unit,* )dtareasic(i,j)
         enddo
      enddo

c  10 format(10f10.4)
      end


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                                  #
! #                                                             #
! #            LAMBERTIAN_FUNCTION                              #
! #            ROSSTHIN_FUNCTION                                #
! #            ROSSTHICK_FUNCTION                               #
! #            LISPARSE_FUNCTION                                #
! #            LIDENSE_FUNCTION                                 #
! #            ROUJEAN_FUNCTION                                 #
! #            HAPKE_FUNCTION                                   #
! #            HAPKE_FUNCTION_OLD                               #
! #            HAPKEKER  (private)                              #
! #            RAHMAN_FUNCTION                                  #
! #            COXMUNK_FUNCTION                                 #
! #            COXMUNK_FUNCTION_DB                              #
! #                                                             #
! #  These two kernels introduced for Version 3.4R              #
! #                                                             #
! #            BREONVEG_FUNCTION                                #
! #            BREONSOIL_FUNCTION                               #
! #                                                             #
! # New Cox-Munk Subroutines in this Module (Version 3.7)       #
! #                                                             #
! #            BRDF_Generalized_glint                           #
! #            BRDF_Water_RefracIndex                           #
! #            BRDF_WhiteCap_Reflectance                        #
! #                                                             #
! ###############################################################

      MODULE brdf_sup_kernels_m

      use brdf_sup_aux_m, only : derfc_e, BRDF_Fresnel_Complex

      PRIVATE HAPKEKER
      PUBLIC :: LAMBERTIAN_FUNCTION,     &
                ROSSTHIN_FUNCTION,       &
                ROSSTHICK_FUNCTION,      &
                LISPARSE_FUNCTION,       &
                LIDENSE_FUNCTION,        &
                ROUJEAN_FUNCTION,        &
                HAPKE_FUNCTION,          &
                HAPKE_FUNCTION_OLD,      &
                RAHMAN_FUNCTION,         &
                COXMUNK_FUNCTION,        &
                COXMUNK_FUNCTION_DB,     &
                BREONVEG_FUNCTION,       &
                BREONSOIL_FUNCTION,      &
                BRDF_Generalized_Glint,  &
                BRDF_Water_RefracIndex,  &
                BRDF_WhiteCap_Reflectance

      CONTAINS

      SUBROUTINE LAMBERTIAN_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, LAMBERTIAN_KERNEL )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ONE

      implicit none

!  Subroutine arguments

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

!  Lambertian kernel

      LAMBERTIAN_KERNEL = ONE

!  Finish

      RETURN
      END SUBROUTINE LAMBERTIAN_FUNCTION

!

      SUBROUTINE ROSSTHIN_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, ROSSTHIN_KERNEL )

!  module, dimensions and numbers

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

      implicit none

!  Subroutine arguments

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

!  Local variables

      REAL(fpk)  :: DS1, DS2, CKSI, SKSI, KSI, FUNC
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise

      ROSSTHIN_KERNEL = ZERO
      XPHI = PIE - PHI
      CKPHI = - CPHI

!  kernel

      DS1 = XI * XJ
      DS2 = SXI * SXJ
      CKSI = DS1 + DS2 * CKPHI
      IF ( CKSI.GT.ONE ) CKSI = ONE
      SKSI = DSQRT(ONE-CKSI*CKSI)
      KSI = DACOS(CKSI)
      FUNC = ((PIO2-KSI)*CKSI + SKSI)/DS1
      ROSSTHIN_KERNEL = FUNC - PIO2

!  Finish

      RETURN
      END SUBROUTINE ROSSTHIN_FUNCTION

!

      SUBROUTINE ROSSTHICK_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, ROSSTHICK_KERNEL )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, ONE, PIE, PIO2, PIO4

      implicit none

!  Subroutine arguments

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

!  Local variables

      REAL(fpk)  :: DS1, DS2, DS3, CKSI, SKSI, KSI, FUNC
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise

      ROSSTHICK_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  kernel

      DS1 = XI * XJ
      DS2 = SXI * SXJ
      DS3 = XI  + XJ
      CKSI = DS1 + DS2 * CKPHI
      IF ( CKSI.GT.ONE ) CKSI = ONE
      SKSI = DSQRT(ONE-CKSI*CKSI)
      KSI = DACOS(CKSI)
      FUNC = ((PIO2-KSI)*CKSI + SKSI)/DS3
      ROSSTHICK_KERNEL = FUNC - PIO4

!  Finish

      RETURN
      END SUBROUTINE ROSSTHICK_FUNCTION

!

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

!  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 )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: LISPARSE_KERNEL

!  local variables

      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)  :: XPHI, CKPHI

!  Initialise
!    -- Return for special case

      LISPARSE_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI
      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

      CKSI = X_INC  * X_REF + SX_INC * SX_REF * CKPHI

!  contributions P and R

!    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

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

!  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

      IF ( COST .GT. ONE ) THEN
        Q = ONE
      ELSE
        T        = DACOS(COST)
        SINT     = DSQRT ( ONE - COST * COST )
        SINTCOST = SINT * COST
        Q = ONE -  ( ( T - SINTCOST ) / PIE )
      ENDIF

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

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

!  Finish

      RETURN
      END SUBROUTINE LISPARSE_FUNCTION

!

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

!  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 )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: LIDENSE_KERNEL

!  local variables

      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)  :: XPHI, CKPHI

!  Initialise
!    -- Return for special case

      LIDENSE_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI
      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

      CKSI = X_INC  * X_REF + SX_INC * SX_REF * CKPHI

!  contributions P and R

!    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

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

!  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

      IF ( COST .GT. ONE ) THEN
        Q = ONE
      ELSE
        T        = DACOS(COST)
        SINT     = DSQRT ( ONE - COST * COST )
        SINTCOST = SINT * COST
        Q = ONE -  ( ( T - SINTCOST ) / PIE )
      ENDIF

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

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

!  Finish

      RETURN
      END SUBROUTINE LIDENSE_FUNCTION

!

      SUBROUTINE ROUJEAN_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, ROUJEAN_KERNEL )

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, TWO, PIE, PIO2, PI4

      implicit none

!  Subroutine arguments

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

!  Local variables

      REAL(fpk)  :: DS1, DS2, DS3, TXJ, TXI, PHIFAC, S1, S2
      REAL(fpk)  :: XPHI_R, CXPHI_R, SXPHI_R, XPHI_C
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise

      ROUJEAN_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  kernel

      XPHI_C = XPHI
      IF ( XPHI .GT. PIE )  XPHI_C = TWO*PIE - XPHI
      IF ( XPHI .LT. ZERO ) XPHI_C = - XPHI

      IF ( SXI .LT. ZERO ) THEN
        XPHI_R  = ( PIE - XPHI_C )
        CXPHI_R = DCOS ( XPHI_R )
        SXPHI_R = DSIN ( XPHI_R )
        TXI =  - ( SXI / XI )
      ELSE
        TXI =   ( SXI / XI )
        XPHI_R  = XPHI_C
        CXPHI_R = DCOS ( XPHI_R )
        SXPHI_R = DSIN ( XPHI_R )
      ENDIF

      TXJ =  ( SXJ / XJ )
      DS1 = TWO * TXJ * TXI
      DS2 = TXJ + TXI
      DS3 = TXJ*TXJ  + TXI*TXI
      PHIFAC = ( ( PIE - XPHI_R ) * CXPHI_R + SXPHI_R ) / PI4
      S1 = PHIFAC * DS1
      S2 = ( DS2 + DSQRT ( DS3 - DS1 * CXPHI_R ) ) / PIE
      ROUJEAN_KERNEL = S1 - S2

!  Finish

      RETURN
      END SUBROUTINE ROUJEAN_FUNCTION

!

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

!  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 )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: HAPKE_KERNEL

!  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

      REAL(fpk)  :: CTHETA, THETA, PHASE
      REAL(fpk)  :: HOTSPOT, B0_EMPIR, HELP_HOT, B_HOT
      REAL(fpk)  :: SSALBEDO, GAMMA, REFLEC, FUNCTION
      REAL(fpk)  :: HELP_J, TERM_J, HELP_I, TERM_I
      REAL(fpk)  :: XPHI, CKPHI

!  Initialise

      HAPKE_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  kernel

!  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
      IF ( CTHETA .GT. ONE ) CTHETA = ONE
      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
      TERM_J   = ( ONE + HELP_J ) / ( ONE + HELP_J * GAMMA )
      HELP_I   = TWO * XI
      TERM_I   = ( ONE + HELP_I ) / ( ONE + HELP_I * GAMMA )

!  Function

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

!  Finish

      RETURN
      END SUBROUTINE HAPKE_FUNCTION

!

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

!  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 )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: RAHMAN_KERNEL

!  Revision. 24 October 2007.
!  --------------------------

!    * Limiting cases and hotspot evaluation.
!    * Revision based on the DISORT_2 code
!    *  In Disort, this kernel is known as the RPV^ BRDF.

!     The RPV reference is:
!       Rahman, Pinty, Verstraete, 1993: Coupled Surface-Atmosphere
!       Reflectance (CSAR) Model. 2. Semiempirical Surface Model Usable
!       With NOAA Advanced Very High Resolution Radiometer Data,
!       J. Geophys. Res., 98, 20791-20801.

!  The hotspot should occur when XI = XJ and PHI = 180.

!  local variables

      REAL(fpk)  :: T_INC, T_REF, DT1, DT2
      REAL(fpk)  :: CXI, DELTA, K1_SQ, FACT
      REAL(fpk)  :: GEOM, PHASE, RFAC, K0, K1, K2
      REAL(fpk)  :: XPHI, CKPHI, HSPOT, UPPER_LIMIT
      REAL(fpk), PARAMETER  :: SMALL = 1.0d-04

!  Initial section
!  ---------------

!  Initialise output

      RAHMAN_KERNEL = ZERO

!  Limiting case, formerly
!      IF ( XI.EQ.ZERO .OR. XJ.EQ.ZERO ) RETURN

!  Limiting case, revised

      IF ( XJ.LT.SMALL ) RETURN

!  Azimuth convettion

      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  parameters

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

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

!  Value of hot spot

      FACT = K0 * ( TWO - K0 )
      FACT = FACT * ( ONE - K1 ) / ( ONE + K1 ) / ( ONE + K1 )
      GEOM = ( TWO * XJ * XJ * XJ ) ** ( K2 - ONE )
      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
      FACT  = ( ONE + K1_SQ + TWO * K1 * CXI ) ** ONEP5
      PHASE = ( ONE - K1_SQ ) / FACT

!  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 )
      RFAC = ( ONE - K0 ) / ( ONE + DELTA )

!  Geom factor and kernel

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

!  Check upper limit not exceeded

      IF ( RAHMAN_KERNEL .GT. UPPER_LIMIT ) THEN
        RAHMAN_KERNEL = UPPER_LIMIT
      ENDIF

!  Finish

      RETURN
      END SUBROUTINE RAHMAN_FUNCTION

!

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

!  module, dimensions and numbers

      USE LIDORT_pars, only : fpk, ZERO, PIE

      implicit none

!  Subroutine arguments

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

!  Hapke Kernel function.
!   From DISORT code; used as validation.

!  local variables

      REAL(fpk)  :: XPHI, CKPHI
      REAL       :: MU, MUP, DUMMY, DPHI

!  Initialise

      HAPKE_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  Kernel

      MUP = SNGL(XJ)
      MU = SNGL(XI)
      DPHI = SNGL(XPHI)
      DUMMY = ZERO
      HAPKE_KERNEL = HAPKEKER( DUMMY, DUMMY, MU, MUP, DPHI )

!  Finish

      RETURN
      END SUBROUTINE HAPKE_FUNCTION_OLD

!

      REAL FUNCTION HAPKEKER( WVNMLO, WVNMHI, MU, MUP, DPHI )

      implicit none

!      Supplies surface bi-directional reflectivity.
!
!      NOTE 1: Bidirectional reflectivity in DISORT is defined
!              by Eq. 39 in STWL.
!      NOTE 2: Both MU and MU0 (cosines of reflection and incidence
!              angles) are positive.
!
!  INPUT:
!
!    WVNMLO : Lower wavenumber (inv cm) of spectral interval
!
!    WVNMHI : Upper wavenumber (inv cm) of spectral interval
!
!    MU     : Cosine of angle of reflection (positive)
!
!    MUP    : Cosine of angle of incidence (positive)
!
!    DPHI   : Difference of azimuth angles of incidence and reflection
!                (radians)
!
!  LOCAL VARIABLES:
!
!    IREF   : bidirectional reflectance options
!             1 - Hapke's BDR model
!
!    B0     : empirical factor to account for the finite size of
!             particles in Hapke's BDR model
!
!    B      : 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
!
!    H0     : H( mu0 ) in Hapke's BDR model
!
!    H      : H( mu ) in Hapke's BDR model
!
!    HH     : angular width parameter of opposition effect in Hapke's
!             BDR model
!
!    P      : scattering phase function in Hapke's BDR model
!
!    THETA  : phase angle (radians); the angle between incidence and
!             reflection directions in Hapke's BDR model
!
!    W      : single scattering albedo in Hapke's BDR model
!
!
!   Called by- DREF, SURFAC
! +-------------------------------------------------------------------+
!     .. Scalar Arguments ..

      REAL      :: DPHI, MU, MUP, WVNMHI, WVNMLO
!     ..
!     .. Local Scalars ..

      INTEGER   :: IREF
      REAL      :: B0, B, CTHETA, GAMMA, H0, H, HH, P, THETA, W
!     ..
!     .. Intrinsic Functions ..

      INTRINSIC COS, SQRT
!     ..

      IREF = 1

      IF ( IREF.EQ.1 ) THEN

!                              ** Hapke's BRDF model (times Pi/Mu0)
!                              ** (Hapke, B., Theory of reflectance
!                              ** and emittance spectroscopy, Cambridge
!                              ** University Press, 1993, Eq. 8.89 on
!                              ** page 233. Parameters are from
!                              ** Fig. 8.15 on page 231, expect for w.)

         CTHETA = MU * MUP + (1.-MU**2)**.5 * (1.-MUP**2)**.5 * COS( DPHI )
         THETA = ACOS( CTHETA )

         P    = 1. + 0.5 * CTHETA

         HH   = 0.06
         B0   = 1.0
         B    = B0 * HH / ( HH + TAN( THETA/2.) )

         W = 0.6
         GAMMA = SQRT( 1. - W )
         H0   = ( 1. + 2.*MUP ) / ( 1. + 2.*MUP * GAMMA )
         H    = ( 1. + 2.*MU ) / ( 1. + 2.*MU * GAMMA )

         hapkeker = W / 4. / (MU+MUP) * ( (1.+B)* P + H0 * H - 1.0 )

      END IF

      RETURN
      END FUNCTION HAPKEKER

!

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

!  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 )
      REAL(fpk), intent(in)  :: XI, SXI, XJ, SXJ, PHI, CPHI, SKPHI
      REAL(fpk), intent(out) :: COXMUNK_KERNEL

!  Critical exponent taken out

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

!  Local variables

      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
      REAL(fpk)  :: S1, S2, S3, XXI, XXJ, T1, T2, DCOT
      REAL(fpk)  :: SHADOWI, SHADOWR, 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
      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 )
      RP = ( H1 - H2 ) / ( H1 + H2 )
      RL = ( Z2 - H2 ) / ( Z2 + H2 )
      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)
      ARGUMENT = TA * TA  / 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

!  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 = XI/DSQRT(ONE-XXI)
      T1   = DEXP(-DCOT*DCOT*S2)
!      T2   = DERFC(DCOT*S3)
      T2   = DERFC_E(DCOT*S3)
      SHADOWI = HALF*(S1*T1/DCOT-T2)

      XXJ  = XJ*XJ
      DCOT = XJ/DSQRT(ONE-XXJ)
      T1   = DEXP(-DCOT*DCOT*S2)
!      T2   = DERFC(DCOT*S3)
      T2   = DERFC_E(DCOT*S3)
      SHADOWR = HALF*(S1*T1/DCOT-T2)

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

!     Finish

      RETURN
      END SUBROUTINE COXMUNK_FUNCTION

!

      SUBROUTINE COXMUNK_FUNCTION_DB  &
      ( MAXPARS, NPARS, 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 )

!  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

!  Subroutine arguments

      INTEGER  , intent(in)  :: MAXPARS, NPARS
      REAL(fpk), intent(in)  :: 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

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

!  help variables

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

!  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)  :: RHOLD       ( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )
      REAL(fpk)  :: R0_MSRS_QUAD( MAX_MSRS_MUQUAD, MAX_MSRS_PHIQUAD )

!  Safety first zeroing

      REFLEC_0 = ZERO
      COXMUNK_KERNEL = ZERO

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

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

!  Higher orders scattering

      REFLEC   = REFLEC_0
      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 &
              ( MAXPARS, NPARS, PARS, &
                XM, SXM, XI, SXI, PHI_SUB2, CPHI_SUB2, SPHI_SUB2, &
                R0_OUT_QUAD(K,N) )
            CALL COXMUNK_FUNCTION &
              ( MAXPARS, NPARS, PARS, &
                XJ, SXJ, XM, SXM, PHI_SUB1, CPHI_SUB1, SPHI_SUB1, &
                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

!  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 &
                       ( MAXPARS, NPARS, PARS, &
                         XMI, SXMI, XMR, SXMR, PHI_SUB1, CPHI_SUB1, SPHI_SUB1, &
                         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

!  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

!  Continuation point for finishing MSR

 67      CONTINUE

!  Add to total

         REFLEC = REFLEC + REFLEC_S

!  End scattering order loop

      ENDDO

!  Compute total

      COXMUNK_KERNEL = REFLEC

!  debug

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

!  Finish

      RETURN
      END SUBROUTINE COXMUNK_FUNCTION_DB

!

      SUBROUTINE BREONVEG_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, BREONVEG_KERNEL )

!  module, dimensions and numbers

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

      implicit none

!  Subroutine arguments

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

!  Local variables

      REAL(fpk)  :: Z, Z1, Z2, Z2_SQ_M1, H1, H2, RP, RL, FP
      REAL(fpk)  :: XPHI, CKPHI, ATTEN, PROJECTIONS
      REAL(fpk)  :: sgamma, cgamma, calpha, calpha_sq, salpha
      REAL(fpk)  :: PLEAF, GS, GV, FP0

!  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 vegetation model (2009)

!  Initialise

      BREONVEG_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  ..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
!   PARS(1) = refractive index squared

      Z2_SQ_M1 = Z2 * Z2 + MINUS_ONE
      H1 = PARS(1) * Z2
      H2 = DSQRT ( PARS(1) + Z2_SQ_M1 )
      RP = ( H1 - H2 ) / ( H1 + H2 )
      RL = ( Z2 - H2 ) / ( Z2 + H2 )
      FP = HALF * ( RP*RP + RL*RL )

!  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    = dsqrt(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.0d0 * calpha_sq * salpha  / pie

! Polarization model for vegetation

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

! BRDF  with attenuation factor

      cgamma = Z2
      sgamma = dsqrt ( one - cgamma * cgamma )
      atten  = one - sgamma
      BREONVEG_KERNEL = Fp0 * atten

!     Finish

      RETURN
      END SUBROUTINE BREONVEG_FUNCTION

!

      SUBROUTINE BREONSOIL_FUNCTION  &
   ( MAXPARS, NPARS, PARS, XJ, SXJ, XI, SXI, PHI, CPHI, SKPHI, BREONSOIL_KERNEL )

!  module, dimensions and numbers

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

      implicit none

!  Subroutine arguments

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

!  Local variables

      REAL(fpk)  :: Z, Z1, Z2, Z2_SQ_M1, H1, H2, RP, RL, FP
      REAL(fpk)  :: XPHI, CKPHI, ATTEN, FP1
      REAL(fpk)  :: sgamma, cgamma, calpha, calpha_sq, salpha

!  F-.M. Breon Soil model (2009)

!  Initialise

      BREONSOIL_KERNEL = ZERO
      XPHI  = PIE - PHI
      CKPHI = - CPHI

!  ..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
!   PARS(1) = refractive index squared

      Z2_SQ_M1 = Z2 * Z2 + MINUS_ONE
      H1 = PARS(1) * Z2
      H2 = DSQRT ( PARS(1) + Z2_SQ_M1 )
      RP = ( H1 - H2 ) / ( H1 + H2 )
      RL = ( Z2 - H2 ) / ( Z2 + H2 )
      FP = HALF * ( RP*RP + RL*RL )

!  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    = dsqrt(one - calpha_sq)

! Polarization model for soil

      Fp1 = 0.25_fpk * FP / xi / xj

! BRDF  with attenuation factor

      cgamma = Z2
      sgamma = dsqrt ( one - cgamma * cgamma )
      atten  = one - sgamma
      BREONSOIL_KERNEL = Fp1 * atten

!     Finish

      RETURN
      END SUBROUTINE BREONSOIL_FUNCTION

!

subroutine BRDF_WhiteCap_Reflectance &
    ( WindSpeed, Wavelength, WC_Reflectance, WC_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)

!   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

!  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, WLP, Ref(39), wspd, wl, Ref_i, Rwc

!  Initialize

   WC_Reflectance = zero
   WC_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
   IF (wspd .le. 9.25) THEN
       Wlb = 0.01*((3.18e-03)*((wspd-3.7)**3.0))
   ELSE IF (wspd .gt. 9.25) THEN
       Wlb = 0.01*((4.82e-04)*((wspd+1.8)**3.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

!  Final values

   WC_Lambertian  = real(Wlb,fpk)
   WC_Reflectance = real(Rwc,fpk)

!  Finish

   return
end subroutine BRDF_WhiteCap_Reflectance

!

subroutine BRDF_Water_RefracIndex &
     ( Wavelength, Salinity, Refrac_R, Refrac_I )

   use lidort_pars, only : fpk, zero

!  THIS IS FORMERLY CALLED "INDWAT"

   implicit none

!  Input/Output
!  ------------

   real(fpk), intent(in)   :: Wavelength
   real(fpk), intent(in)   :: Salinity
   real(fpk), intent(out)  :: Refrac_R
   real(fpk), intent(out)  :: Refrac_I

!  Local (Real-valued quantities)
!  -----

!subroutine indwat(wl,xsal,nr,ni)
!
! input parameters:  wl=wavelength (in micrometers)
!                    xsal=salinity (in ppt), if xsal<0 then 34.3ppt by default
! output parameters: nr=index of refraction of sea water
!                    ni=extinction coefficient of sea water

       real twl(62),tnr(62),tni(62)
       real nr,ni,wl,xwl,yr,yi,nrc,nic,xsal
       integer i

! Indices of refraction for pure water from Hale and Querry,
! Applied Optics, March 1973, Vol. 12,  No. 3, pp. 555-563
       data twl/&
        0.250,0.275,0.300,0.325,0.345,0.375,0.400,0.425,0.445,0.475,&
        0.500,0.525,0.550,0.575,0.600,0.625,0.650,0.675,0.700,0.725,&
        0.750,0.775,0.800,0.825,0.850,0.875,0.900,0.925,0.950,0.975,&
        1.000,1.200,1.400,1.600,1.800,2.000,2.200,2.400,2.600,2.650,&
        2.700,2.750,2.800,2.850,2.900,2.950,3.000,3.050,3.100,3.150,&
        3.200,3.250,3.300,3.350,3.400,3.450,3.500,3.600,3.700,3.800,&
        3.900,4.000/
        data tnr/&
        1.362,1.354,1.349,1.346,1.343,1.341,1.339,1.338,1.337,1.336,&
        1.335,1.334,1.333,1.333,1.332,1.332,1.331,1.331,1.331,1.330,&
        1.330,1.330,1.329,1.329,1.329,1.328,1.328,1.328,1.327,1.327,&
        1.327,1.324,1.321,1.317,1.312,1.306,1.296,1.279,1.242,1.219,&
        1.188,1.157,1.142,1.149,1.201,1.292,1.371,1.426,1.467,1.483,&
        1.478,1.467,1.450,1.432,1.420,1.410,1.400,1.385,1.374,1.364,&
        1.357,1.351/
        data tni/&
        3.35E-08,2.35E-08,1.60E-08,1.08E-08,6.50E-09,&
        3.50E-09,1.86E-09,1.30E-09,1.02E-09,9.35E-10,&
        1.00E-09,1.32E-09,1.96E-09,3.60E-09,1.09E-08,&
        1.39E-08,1.64E-08,2.23E-08,3.35E-08,9.15E-08,&
        1.56E-07,1.48E-07,1.25E-07,1.82E-07,2.93E-07,&
        3.91E-07,4.86E-07,1.06E-06,2.93E-06,3.48E-06,&
        2.89E-06,9.89E-06,1.38E-04,8.55E-05,1.15E-04,&
        1.10E-03,2.89E-04,9.56E-04,3.17E-03,6.70E-03,&
        1.90E-02,5.90E-02,1.15E-01,1.85E-01,2.68E-01,&
        2.98E-01,2.72E-01,2.40E-01,1.92E-01,1.35E-01,&
        9.24E-02,6.10E-02,3.68E-02,2.61E-02,1.95E-02,&
        1.32E-02,9.40E-03,5.15E-03,3.60E-03,3.40E-03,&
        3.80E-03,4.60E-03/

!  Assign input

      wl   = real(WAVELENGTH)
      xsal = real(SALINITY)
      Refrac_R = zero
      Refrac_I = zero

!  Find wavelength point for interpolation

        i=2
 10     if (wl.lt.twl(i)) goto 20
        if (i.lt.62) then
           i=i+1
           goto 10
         endif

!  Interpolate

 20     xwl=twl(i)-twl(i-1)
        yr=tnr(i)-tnr(i-1)
        yi=tni(i)-tni(i-1)
        nr=tnr(i-1)+(wl-twl(i-1))*yr/xwl
        ni=tni(i-1)+(wl-twl(i-1))*yi/xwl
!
! Correction to be applied to the index of refraction and to the extinction
! coefficients of the pure water to obtain the ocean water one (see for
! example Friedman). By default, a typical sea water is assumed
! (Salinity=34.3ppt, Chlorinity=19ppt) as reported by Sverdrup.
! In that case there is no correction for the extinction coefficient between
! 0.25 and 4 microns. For the index of refraction, a correction of +0.006
! has to be applied (McLellan). For a chlorinity of 19.0ppt the correction
! is a linear function of the salt concentration. Then, in 6S users are able
! to enter the salt concentration (in ppt).
! REFERENCES:
! Friedman D., Applied Optics, 1969, Vol.8, No.10, pp.2073-2078.
! McLellan H.J., Elements of physical Oceanography, Pergamon Press, Inc.,
!        New-York, 1965, p 129.
! Sverdrup H.V. et al., The Oceans (Prentice-Hall, Inc., Englewood Cliffs,
!        N.J., 1942, p 173.

        nrc=0.006
        nic=0.000
        nr=nr+nrc*(xsal/34.3)
        ni=ni+nic*(xsal/34.3)

!  Assign output

    REFRAC_R = real(nr,fpk)
    REFRAC_I = real(ni,fpk)

    return
end subroutine BRDF_Water_RefracIndex

SUBROUTINE BRDF_Generalized_Glint &
         ( 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, SUNGLINT_REFLEC )

     use lidort_pars, only : fpk, zero, one, half, quarter, two, three, four, 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

!  Cox-Munk Coefficients. Intent(inout).

      real(fpk), intent(inout) :: SUNGLINT_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

!  Initialise output

      SUNGLINT_REFLEC = ZERO

!  COmpute coefficients, according to 6S formulation

      IF ( DO_COEFFS ) THEN
         SUNGLINT_COEFFS = zero
         IF ( DO_ISOTROPIC ) THEN
            SUNGLINT_COEFFS(1) = 0.003_fpk + 0.00512_fpk * WINDSPEED
         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
         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))
         WSigU = ONE / Sqrt(SUNGLINT_COEFFS(2))
         VAR   = WSigC * WSigU * HALF

!  angles

         XE = (  CKPHI_W * ZX + SKPHI_W * ZY ) * WSigC ; XE_sq = XE * XE ; XE_sq_1 = xe_sq - one
         XN = ( -SKPHI_W * ZX + CKPHI_W * ZY ) * WSigU ; XN_sq = XN * XN ; XN_sq_1 = xn_sq - one

!  GC Coefficient

         Coeff = ONE - SUNGLINT_COEFFS(3) *      XE_sq_1      * XN * half &
                     - SUNGLINT_COEFFS(4) * ( XN_sq - three ) * XN / six  &
                     + SUNGLINT_COEFFS(5) * ( XE_sq * XE_sq - six * XE_sq + three ) / twentyfour &
                     + SUNGLINT_COEFFS(6) * ( XN_sq * XN_sq - six * XN_sq + three ) / twentyfour &
                     + SUNGLINT_COEFFS(7) * XE_sq_1 * XN_sq_1 / four

!  Probability and finish

         ARGUMENT = ( XE_sq  + XN_sq ) * HALF
         IF ( ARGUMENT .LT. CRITEXP ) THEN
            PROB = COEFF * EXP ( - ARGUMENT ) * VAR
            FAC2 = QUARTER / XI / XJ / ( COSTILT ** FOUR )
            SUNGLINT_REFLEC = XMP * PROB * FAC2
         ENDIF

      ENDIF

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

      IF ( DO_ISOTROPIC ) THEN

!  Compute Probability and finish

         VAR   = SUNGLINT_COEFFS(1)
         ARGUMENT = TANTILT_SQ / VAR
         IF ( ARGUMENT .LT. CRITEXP ) THEN
            PROB = EXP ( - ARGUMENT ) / VAR
            FAC2 = QUARTER / XI / XJ / ( COSTILT ** FOUR )
            SUNGLINT_REFLEC = XMP * PROB * 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

!  End module

      END MODULE brdf_sup_kernels_m

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps