text:
! MD5 of template: 8ca8420e02fee968b6f8d829e4f8b5c7
! Array related routines (Integration, Interpolation, etc.)
!
! ------------------------------------------------------------------------------
! Copyright (c) 2009-13, Thomas P. Robitaille
!
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
!  * Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
!
!  * Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ------------------------------------------------------------------------------
!
! This module contains subroutines to make it easy to read in
! configuration files. The expected format of a .conf file is
! an ASCII file with lines containing lines such as:
!
! parameter = value
!
! Note that character variables should be enclosed in inverted commas
! especially when giving paths to files. Commented lines (starting with
! a #) and blank lines are allowed. Comments directly following a
! parameter statement on the same line are NOT allowed:
!
! parameter = value  # this is an INVALID comment
!
! The subroutines are:
!
! load_config_file(filename) - load a config file into RAM, ignoring
!                              comments and blank lines
!
! read_config(parameter,value) - read in the value for 'parameter'
!
! read_config is in fact an interface to four different subroutines
! depending on the type of the 'value' variable. This variable can
! be real(sp), real(dp), integer, or character(len=*)

module lib_conf

  implicit none
  save

  private
  public :: load_config_file

  integer,parameter :: idp = selected_int_kind(13)
  integer,parameter :: sp = selected_real_kind(p=6,r=37)
  integer,parameter :: dp = selected_real_kind(p=15,r=307)

  integer,parameter :: config_file_len=1000
  integer,parameter :: config_line_len=100
  character(len=6),parameter :: config_line_fmt='(A100)'

  integer,private :: n_lines
  ! number of lines is the currently loaded config file

  character(len=config_line_len),dimension(config_file_len),private :: line
  ! the lines of the config file

  public :: read_config
  interface read_config
     module procedure config_real8,config_real4,config_char,config_int4,config_int8,config_logical
  end interface read_config

  interface read_from_string
     module procedure read_int4_from_string
     module procedure read_int8_from_string
     module procedure read_real4_from_string
     module procedure read_real8_from_string
     module procedure read_char_from_string
     module procedure read_logical_from_string
  end interface read_from_string

contains

  subroutine read_int4_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    integer,intent(out) :: value
    read(string,*) value
  end subroutine read_int4_from_string

  subroutine read_int8_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    integer(idp),intent(out) :: value
    read(string,*) value
  end subroutine read_int8_from_string

  subroutine read_real4_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    real(sp),intent(out) :: value
    read(string,*) value
  end subroutine read_real4_from_string

  subroutine read_real8_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    real(dp),intent(out) :: value
    read(string,*) value
  end subroutine read_real8_from_string

  subroutine read_char_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    character(len=*),intent(out) :: value
    value = trim(adjustl(string))
  end subroutine read_char_from_string

  subroutine read_logical_from_string(string,value)
    implicit none
    character(len=*),intent(in) :: string
    logical,intent(out) :: value
    select case(trim(adjustl(string)))
    case('Y','y','yes','YES','Yes','true',"'YES'")
       value=.true.
    case('N','n','no','NO','No','false',"'NO'")
       value=.false.
    case default
       print *,"unknown logical"//trim(string)
    end select
  end subroutine read_logical_from_string

  subroutine load_config_file(filename)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: filename
    ! the name of the .conf file to read in

    ! --- Local variables --- !

    character(len=config_line_len) :: line_temp

    integer :: ioerr
    ! used for I/O errors

    open(unit=80, file=filename, status='old')

    n_lines = 0
    line = ''

    do

       ! --- Read in the next line --- !

       read(80,config_line_fmt,iostat=ioerr) line_temp
       if(ioerr.ne.0) exit

       ! --- Remove preceding whitespace --- !

       line_temp = adjustl(line_temp)

       ! --- Ignore lines with nothing and lines with comments --- !

       if(line_temp(1:1).ne.'#'.and.trim(line_temp).ne.'') then

          n_lines = n_lines + 1
          line(n_lines) = line_temp

       end if

    end do

    close(unit=80)

  end subroutine load_config_file

  subroutine config_logical(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    logical,intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_logical

  subroutine config_real8(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    real(dp),intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_real8

  subroutine config_real4(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    real(sp),intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_real4

  subroutine config_char(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    character(len=*),intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_char

  subroutine config_int8(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    integer(idp),intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_int8

  subroutine config_int4(par_name,value,element)

    implicit none

    ! --- Input --- !

    character(len=*),intent(in) :: par_name
    ! parameter to search for

    integer,optional,intent(in) :: element
    ! array element (if needed)

    ! --- Output --- !

    integer,intent(out) :: value
    ! the value that was read in

    character(len=100) :: c_element,par_name_new,name_check
    ! temporary parameter name (e.g. with array index)

    integer :: i,pos
    ! loop and position variables

    logical :: found
    ! whether the parameter was found

    found = .false.

    if(present(element)) then
       write(c_element,'(I0)') element
       par_name_new = trim(par_name)//'('//trim(c_element)//')'
    else
       par_name_new = par_name
    end if

    do i=1,n_lines

       if(index(line(i),trim(par_name_new)).gt.0) then

          pos = index(line(i),'=')
          name_check = adjustl(line(i)(1:pos-1))
          if(trim(name_check) == trim(par_name_new)) then
             call read_from_string(line(i)(pos+1:),value)
             found = .true.
             exit
          end if

       end if

    end do

    if(.not.found) then
       print *, "Parameter not found : "//trim(par_name_new)
       stop
    end if

  end subroutine config_int4

end module lib_conf

max_stars_repo_name: jfitz/code-stat

max_issues_repo_name: jfitz/code-stat

max_forks_repo_name: jfitz/code-stat


text:
program sample_coupler
  use mpi
  use ctca
  implicit none

  integer(kind=4) :: dataintnum, myrank, nprocs, ierr, areaid, progid, fromrank
  integer(kind=4),dimension(4) :: reqinfo
  integer(kind=4),dimension(2) :: dataint
  real(kind=8),dimension(6,400) :: datareal8

  call CTCAC_init()
  print *, "coupler init done"

  call MPI_Comm_size(CTCA_subcomm, nprocs, ierr)
  call MPI_Comm_rank(CTCA_subcomm, myrank, ierr)

  call CTCAC_regarea_real8(areaid)

  do while (.true.)
     call CTCAC_pollreq(reqinfo, fromrank, dataint, 2)
  print *, "coupler poll req done"

     if (CTCAC_isfin()) then
  print *, "coupler got fin"
        exit
     end if

     call CTCAC_readarea_real8(areaid, fromrank, (dataint(2)-1)*6*400, 6*400, datareal8)
  print *, "coupler read data"

     if (fromrank >= 0) then
  print *, "coupler got req"
        ! decide progid
        progid = dataint(1)

        ! enqueue request with data
        call CTCAC_enqreq_withreal8(reqinfo, progid, dataint, 2, datareal8, 6*400)
  print *, "coupler enqreq done"

     end if
  end do

  call CTCAC_finalize()

end program

max_stars_repo_name: tnanri/cotocoa

max_issues_repo_name: tnanri/cotocoa

max_forks_repo_name: tnanri/cotocoa


text:
c-----------------------------------------------------------------------
c  nek5000 user-file template
c
c  user specified routines:
c     - uservp  : variable properties
c     - userf   : local acceleration term for fluid
c     - userq   : local source term for scalars
c     - userbc  : boundary conditions
c     - useric  : initial conditions
c     - userchk : general purpose routine for checking errors etc.
c     - userqtl : thermal divergence for lowMach number flows
c     - usrdat  : modify element vertices
c     - usrdat2 : modify mesh coordinates
c     - usrdat3 : general purpose routine for initialization
c
c-----------------------------------------------------------------------
      subroutine uservp(ix,iy,iz,eg) ! set variable properties

c      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
c     e = gllel(eg)

      udiff  = 0.0
      utrans = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine userf(ix,iy,iz,eg) ! set acceleration term
c
c     Note: this is an acceleration term, NOT a force!
c     Thus, ffx will subsequently be multiplied by rho(x,t).
c
c      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
c     e = gllel(eg)

      Fr2 = 0.01
      U = 0.125
      f = -1.
      ffx = f*uy
      ffy = -f*ux + f*U
      ffz = temp/Fr2

      return
      end
c-----------------------------------------------------------------------
      subroutine userq(ix,iy,iz,eg) ! set source term

c      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e
c     e = gllel(eg)

c If no stratification restoration
      qvol   = 0.0

c If include stratification restoration
c      qvol = -(1.0/(6.0*pi))*(temp-z)

      return
      end
c-----------------------------------------------------------------------
      subroutine userbc(ix,iy,iz,iside,eg) ! set up boundary conditions
c
c     NOTE ::: This subroutine MAY NOT be called by every process
c
c      implicit none

      integer ix,iy,iz,iside,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

c      if (cbc(iside,gllel(eg),ifield).eq.'v01')

      ux   = 0.0
      uy   = 0.0
      uz   = 0.0
      temp = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine useric(ix,iy,iz,eg) ! set up initial conditions

c      implicit none

      integer ix,iy,iz,eg

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      U = 0.125

      ux   = U
      uy   = 0.0
      uz   = 0.0
      temp = z

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk()

c      implicit none

      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
      include 'ZPER'

      return
      end
c-----------------------------------------------------------------------
      subroutine userqtl ! Set thermal divergence

      call userqtl_scig

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat()   ! This routine to modify element vertices

c      implicit none

      include 'SIZE'
      include 'TOTAL'

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2()  ! This routine to modify mesh coordinates

c      implicit none

      include 'SIZE'
      include 'TOTAL'

      ntot = nx1*ny1*nz1*nelt

      J = 0.6
      chi = 0.16
      U = 0.125
      N = 10.
      f = 1.

      k = N*chi/U
      h = J*U/N

      do i=1,ntot
         xx = xm1(i,1,1,1)
         argx = k*xm1(i,1,1,1)/2.
         argsin = sin(argx)
         zm1(i,1,1,1) = zm1(i,1,1,1) + (1.0-zm1(i,1,1,1))*h*(argsin**2)
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3()

c      implicit none

      include 'SIZE'
      include 'TOTAL'

      return
      end

c automatically added by makenek
      subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
      integer*8 glo_num(1)

      return
      end

ext: f


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! Tests F2003 defined I/O (unformatted write)
! This example is based on example from the Fortran 2003 Handbook (page 336)

module rational_stuff
  type rational
    integer n, d
    real :: rat_rslt=0.0
  end type

  interface write (unformatted)
    module procedure write_rational_value
  end interface write (unformatted)
  interface read (unformatted)
    module procedure read_rational_value
  end interface read (unformatted)

 contains
  subroutine write_rational_value(dtv, unit, iostat, iomsg)

  class(rational), intent(inout) :: dtv
  integer, intent(in) :: unit
  integer, intent(out) :: iostat
  character(len=*), intent(inout) :: iomsg
  real rat

     if (dtv%d .ne. 0) then
       dtv%rat_rslt = real(dtv%n)/real(dtv%d)
       write(unit), dtv%n, dtv%d, dtv%rat_rslt
     endif
  end subroutine

  subroutine read_rational_value(dtv, unit, iostat, iomsg)

  class(rational), intent(inout) :: dtv
  integer, intent(in) :: unit
  integer, intent(out) :: iostat
  character(len=*), intent(inout) :: iomsg

  end subroutine
end module

  use rational_stuff
  type(rational) x, y
  logical rslt(3), expect(3)
  x = rational(2,3)
  open(10, file='io13.output', status='replace', form='unformatted')
  write(10) , x
  close(10)
  open(10, file='io13.output', status='old', form='unformatted')
  read(10), y%n, y%d, y%rat_rslt
  close(10)

  !print *, x%n, x%d, x%rat_rslt
  !print *, y%n, y%d, y%rat_rslt

  expect = .true.
  rslt(1) = x%n .eq. y%n
  rslt(2) = x%d .eq. y%d
  rslt(3) = x%rat_rslt .eq. y%rat_rslt

  call check(rslt, expect, 3)
  end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
!!!!!
! This subroutine performs a calculation of the incomplete Gamma
! function required for calculating the transmission of isotropic
! radiation. For more information on the derivation, see Pinty et al.
! (2006) equations (16) and (17) and surrounding text.
!
!  Pinty, B., Lavergne, T., Dickinson, R.E., Widlowski, J.-L., Gobron,
!  N., Verstraete, M.M., 2006. Simplifying the interaction of land
!  surfaces with radiation for relating remote sensing products to
!  climate models. J. Geophys. Res. 111, D02116.
!  doi:10.1029/2005JD005952
!!!!!

module mod_biggamma
contains
function BigGamma(tau)
    implicit none
    real*8, intent(in) :: tau
    real*8 :: BigGamma

    integer :: order, i

    order = 20
    BigGamma = 0.0d0

    do i=order, 1, -1
        BigGamma = i / (1.0d0 + i / (tau + BigGamma))
    end do

    BigGamma = 1.0d0 / (tau + BigGamma)
end function BigGamma
end module mod_biggamma

max_stars_repo_name: ankurdesai/pecan

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_name: ankurdesai/pecan

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_name: ankurdesai/pecan

max_forks_repo_licenses: NCSA Unlicense


SUBROUTINE READ_CONFIG
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!  Reads in a restart file in which the fluid is a
!!!!  circular droplet
!!!!  K.Travis 24/03/12 Univ. of Sheffield
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
USE header_file
IMPLICIT NONE

REAL(KIND = DOUBLE) :: DUMMY,DET
INTEGER :: I

OPEN(4,FILE = 'restart.txt')

DO I = 1,N
   READ(4,*) X(I),Y(I),Z(I)
   READ(4,*) PX(I),PY(I),PZ(I)
   READ(4,*) MASS(I)
END DO
READ(4,*) EFC(1,1), EFC(1,2),EFC(1,3)
READ(4,*) EFC(2,1), EFC(2,2),EFC(2,3)
READ(4,*) EFC(3,1), EFC(3,2),EFC(3,3)
READ(4,*) ECF(1,1), ECF(1,2),ECF(1,3)
READ(4,*) ECF(2,1), ECF(2,2),ECF(2,3)
READ(4,*) ECF(3,1), ECF(3,2),ECF(3,3)
READ(4,*) ALEN,BLEN,CLEN
READ(4,*) ALPHA,BETA,GAMMA
READ(4,*) DELTA

VOLUME = EFC(1,1)*EFC(2,2)*EFC(3,3)

XYZ_FILE_STAT = 'UNKNOWN'

CLOSE(4)

END SUBROUTINE


text:
C***********************************************************************
      SUBROUTINE SPLINE (N,X,Y,FDP)
C***********************************************************************
C-----THIS SUBROUTINE COMPUTES THE SECOND DERIVATIVES NEEDED
C-----IN CUBIC SPLINE INTERPOLATION.  THE INPUT DATA ARE:
C-----N = NUMBER OF DATA POINTS
C-----X = ARRAY CONTAINING THE VALUES OF THE INDEPENDENT VARIABLE
C-----    (ASSUMED TO BE IN ASCENDING ORDER)
C-----Y = ARRAY CONTAINING THE VALUES OF THE FUNCTION AT THE
C-----    DATA POINTS GIVEN IN THE X ARRAY
C-----THE OUTPUT IS THE ARRAY FDP WHICH CONTAINS THE SECOND
C-----DERIVATIVES OF THE INTERPOLATING CUBIC SPLINE.
      DIMENSION X(N),Y(N),A(N),B(N),C(N),R(N),FDP(N)
C-----COMPUTE THE COEFFICIENTS AND THE RHS OF THE EQUATIONS.
C-----THIS ROUTINE USES THE CANTILEVER CONDITION.  THE PARAMETER
C-----ALAMDA (LAMBDA) IS SET TO 1. BUT THIS CAN BE USER-MODIFIED.
C-----A,B,C ARE THE THREE DIAGONALS OF THE TRIDIAGONAL SYSTEM;
C-----R IS THE RIGHT HAND SIDE.  THESE ARE NOW ASSEMBLED.
      ALAMDA = 1.
      NM2 = N - 2
      NM1 = N - 1
      C(1) = X(2) - X(1)
      DO 1 I=2,NM1
      C(I) = X(I+1) - X(I)
      A(I) = C(I-1)
      B(I) = 2.*(A(I) + C(I))
      R(I) = 6.*((Y(I+1) - Y(I))/C(I) - (Y(I) - Y(I-1))/C(I-1))
    1 CONTINUE
      B(2) = B(2) + ALAMDA * C(1)
      B(NM1) = B(NM1) + ALAMDA * C(NM1)
C-----AT THIS POINT WE COULD CALL A TRIDIAGONAL SOLVER SUBROUTINE
C-----BUT THE NOTATION IS CLUMSY SO WE WILL SOLVE DIRECTLY.  THE
C-----NEXT SECTION SOLVES THE SYSTEM WE HAVE JUST SET UP.
      DO 2 I=3,NM1
      T = A(I)/B(I-1)
      B(I) = B(I) - T * C(I-1)
      R(I) = R(I) - T * R(I-1)
    2 CONTINUE
      FDP(NM1) = R(NM1)/B(NM1)
      DO 3 I=2,NM2
      NMI = N - I
      FDP(NMI) = (R(NMI) - C(NMI)*FDP(NMI+1))/B(NMI)
    3 CONTINUE
      FDP(1) = ALAMDA * FDP(2)
      FDP(N) = ALAMDA * FDP(NM1)
C-----WE NOW HAVE THE DESIRED DERIVATIVES SO WE RETURN TO THE
C-----MAIN PROGRAM.
      RETURN
      END
C***********************************************************************
      SUBROUTINE SPEVAL (N,X,Y,FDP,XX,F)
C***********************************************************************
C-----THIS SUBROUTINE EVALUATES THE CUBIC SPLINE GIVEN
C-----THE DERIVATIVE COMPUTED BY SUBROUTINE SPLINE.
C-----THE INPUT PARAMETERS N,X,Y,FDP HAVE THE SAME
C-----MEANING AS IN SPLINE.
C-----XX = VALUE OF INDEPENDENT VARIABLE FOR WHICH
C-----     AN INTERPOLATED VALUE IS REQUESTED
C-----F =  THE INTERPOLATED RESULT
      DIMENSION X(N),Y(N),FDP(N)
C-----THE FIRST JOB IS TO FIND THE PROPER INTERVAL.
      NM1 = N - 1
      DO 1 I=1,NM1
      IF (XX.LE.X(I+1)) GO TO 10
    1 CONTINUE
C-----NOW EVALUATE THE CUBIC
   10 DXM = XX - X(I)
      DXP = X(I+1) - XX
      DEL = X(I+1) - X(I)
      F = FDP(I)*DXP*(DXP*DXP/DEL - DEL)/6.
     1   +FDP(I+1)*DXM*(DXM*DXM/DEL - DEL)/6.
     2   +Y(I)*DXP/DEL + Y(I+1)*DXM/DEL
      RETURN
      END


text:
! RUN: bbc -emit-fir %s -o - | FileCheck %s

module max_test
    contains
    ! CHECK-LABEL: func @_QMmax_testPdynamic_optional(
    ! CHECK-SAME:  %[[VAL_0:.*]]: !fir.box<!fir.array<?xi32>> {fir.bindc_name = "a"},
    ! CHECK-SAME:  %[[VAL_1:.*]]: !fir.box<!fir.array<?xi32>> {fir.bindc_name = "b"},
    ! CHECK-SAME:  %[[VAL_2:.*]]: !fir.box<!fir.array<?xi32>> {fir.bindc_name = "c", fir.optional}) {
    subroutine dynamic_optional(a, b, c)
      integer :: a(:), b(:)
      integer, optional :: c(:)
    ! CHECK:  %[[VAL_10:.*]] = fir.array_load %[[VAL_0]] : (!fir.box<!fir.array<?xi32>>) -> !fir.array<?xi32>
    ! CHECK:  %[[VAL_11:.*]] = fir.array_load %[[VAL_1]] : (!fir.box<!fir.array<?xi32>>) -> !fir.array<?xi32>
    ! CHECK:  %[[VAL_12:.*]] = fir.is_present %[[VAL_2]] : (!fir.box<!fir.array<?xi32>>) -> i1
    ! CHECK:  %[[VAL_17:.*]] = arith.select %[[VAL_12]], %[[VAL_2]], %{{.*}} : !fir.box<!fir.array<?xi32>>
    ! CHECK:  %[[VAL_18:.*]] = fir.array_load %[[VAL_17]] {fir.optional} : (!fir.box<!fir.array<?xi32>>) -> !fir.array<?xi32>
    ! CHECK:  fir.do_loop %[[VAL_26:.*]] = %{{.*}} to %{{.*}} step %{{.*}} unordered iter_args(%[[VAL_27:.*]] = %{{.*}}) -> (!fir.array<?xi32>) {
    ! CHECK:    %[[VAL_28:.*]] = fir.array_fetch %[[VAL_10]], %[[VAL_26]] : (!fir.array<?xi32>, index) -> i32
    ! CHECK:    %[[VAL_29:.*]] = fir.array_fetch %[[VAL_11]], %[[VAL_26]] : (!fir.array<?xi32>, index) -> i32
    ! CHECK:    %[[VAL_30:.*]] = arith.cmpi sgt, %[[VAL_28]], %[[VAL_29]] : i32
    ! CHECK:    %[[VAL_31:.*]] = arith.select %[[VAL_30]], %[[VAL_28]], %[[VAL_29]] : i32
    ! CHECK:    %[[VAL_32:.*]] = fir.if %[[VAL_12]] -> (i32) {
    ! CHECK:      %[[VAL_33:.*]] = fir.array_fetch %[[VAL_18]], %[[VAL_26]] : (!fir.array<?xi32>, index) -> i32
    ! CHECK:      %[[VAL_34:.*]] = arith.cmpi sgt, %[[VAL_31]], %[[VAL_33]] : i32
    ! CHECK:      %[[VAL_35:.*]] = arith.select %[[VAL_34]], %[[VAL_31]], %[[VAL_33]] : i32
    ! CHECK:      fir.result %[[VAL_35]] : i32
    ! CHECK:    } else {
    ! CHECK:      fir.result %[[VAL_31]] : i32
    ! CHECK:    }
    ! CHECK:    %[[VAL_36:.*]] = fir.array_update %[[VAL_27]], %[[VAL_32]], %[[VAL_26]] : (!fir.array<?xi32>, i32, index) -> !fir.array<?xi32>
    ! CHECK:    fir.result %[[VAL_36]] : !fir.array<?xi32>
    ! CHECK:  }
      print *, max(a, b, c)
    end subroutine

    ! CHECK-LABEL: func @_QMmax_testPdynamic_optional_array_expr_scalar_optional(
    ! CHECK-SAME:  %[[VAL_0:.*]]: !fir.box<!fir.array<?xi32>> {fir.bindc_name = "a"},
    ! CHECK-SAME:  %[[VAL_1:.*]]: !fir.box<!fir.array<?xi32>> {fir.bindc_name = "b"},
    ! CHECK-SAME:  %[[VAL_2:.*]]: !fir.ref<i32> {fir.bindc_name = "c", fir.optional}) {
    subroutine dynamic_optional_array_expr_scalar_optional(a, b, c)
      integer :: a(:), b(:)
      integer, optional :: c
      print *, max(a, b, c)
    ! CHECK:  %[[VAL_10:.*]] = fir.array_load %[[VAL_0]] : (!fir.box<!fir.array<?xi32>>) -> !fir.array<?xi32>
    ! CHECK:  %[[VAL_11:.*]] = fir.array_load %[[VAL_1]] : (!fir.box<!fir.array<?xi32>>) -> !fir.array<?xi32>
    ! CHECK:  %[[VAL_12:.*]] = fir.is_present %[[VAL_2]] : (!fir.ref<i32>) -> i1
    ! CHECK:  fir.do_loop %[[VAL_20:.*]] = %{{.*}} to %{{.*}} step %{{.*}} unordered iter_args(%[[VAL_21:.*]] = %{{.*}}) -> (!fir.array<?xi32>) {
    ! CHECK:    %[[VAL_22:.*]] = fir.array_fetch %[[VAL_10]], %[[VAL_20]] : (!fir.array<?xi32>, index) -> i32
    ! CHECK:    %[[VAL_23:.*]] = fir.array_fetch %[[VAL_11]], %[[VAL_20]] : (!fir.array<?xi32>, index) -> i32
    ! CHECK:    %[[VAL_24:.*]] = arith.cmpi sgt, %[[VAL_22]], %[[VAL_23]] : i32
    ! CHECK:    %[[VAL_25:.*]] = arith.select %[[VAL_24]], %[[VAL_22]], %[[VAL_23]] : i32
    ! CHECK:    %[[VAL_26:.*]] = fir.if %[[VAL_12]] -> (i32) {
    ! CHECK:      %[[VAL_27:.*]] = fir.load %[[VAL_2]] : !fir.ref<i32>
    ! CHECK:      %[[VAL_28:.*]] = arith.cmpi sgt, %[[VAL_25]], %[[VAL_27]] : i32
    ! CHECK:      %[[VAL_29:.*]] = arith.select %[[VAL_28]], %[[VAL_25]], %[[VAL_27]] : i32
    ! CHECK:      fir.result %[[VAL_29]] : i32
    ! CHECK:    } else {
    ! CHECK:      fir.result %[[VAL_25]] : i32
    ! CHECK:    }
    ! CHECK:    %[[VAL_30:.*]] = fir.array_update %[[VAL_21]], %[[VAL_26]], %[[VAL_20]] : (!fir.array<?xi32>, i32, index) -> !fir.array<?xi32>
    ! CHECK:    fir.result %[[VAL_30]] : !fir.array<?xi32>
    ! CHECK:  }
    end subroutine

    ! CHECK-LABEL: func @_QMmax_testPdynamic_optional_scalar(
    ! CHECK-SAME:  %[[VAL_0:.*]]: !fir.ref<i32> {fir.bindc_name = "a"},
    ! CHECK-SAME:  %[[VAL_1:.*]]: !fir.ref<i32> {fir.bindc_name = "b"},
    ! CHECK-SAME:  %[[VAL_2:.*]]: !fir.ref<i32> {fir.bindc_name = "c", fir.optional}) {
    subroutine dynamic_optional_scalar(a, b, c)
      integer :: a, b
      integer, optional :: c
      print *, max(a, b, c)
    ! CHECK:  %[[VAL_8:.*]] = fir.load %[[VAL_0]] : !fir.ref<i32>
    ! CHECK:  %[[VAL_9:.*]] = fir.load %[[VAL_1]] : !fir.ref<i32>
    ! CHECK:  %[[VAL_10:.*]] = fir.is_present %[[VAL_2]] : (!fir.ref<i32>) -> i1
    ! CHECK:  %[[VAL_11:.*]] = arith.cmpi sgt, %[[VAL_8]], %[[VAL_9]] : i32
    ! CHECK:  %[[VAL_12:.*]] = arith.select %[[VAL_11]], %[[VAL_8]], %[[VAL_9]] : i32
    ! CHECK:  %[[VAL_13:.*]] = fir.if %[[VAL_10]] -> (i32) {
    ! CHECK:    %[[VAL_14:.*]] = fir.load %[[VAL_2]] : !fir.ref<i32>
    ! CHECK:    %[[VAL_15:.*]] = arith.cmpi sgt, %[[VAL_12]], %[[VAL_14]] : i32
    ! CHECK:    %[[VAL_16:.*]] = arith.select %[[VAL_15]], %[[VAL_12]], %[[VAL_14]] : i32
    ! CHECK:    fir.result %[[VAL_16]] : i32
    ! CHECK:  } else {
    ! CHECK:    fir.result %[[VAL_12]] : i32
    ! CHECK:  }
    ! CHECK:  fir.call @_FortranAioOutputInteger32(%{{.*}}, %[[VAL_13]]) : (!fir.ref<i8>, i32) -> i1
    end subroutine

    ! CHECK-LABEL: func @_QMmax_testPdynamic_optional_weird(
    ! CHECK-SAME:  %[[VAL_0:.*]]: !fir.ref<i32> {fir.bindc_name = "a"},
    ! CHECK-SAME:  %[[VAL_1:.*]]: !fir.ref<i32> {fir.bindc_name = "b"},
    ! CHECK-SAME:  %[[VAL_2:.*]]: !fir.ref<i32> {fir.bindc_name = "c", fir.optional},
    ! CHECK-SAME:  %[[VAL_3:.*]]: !fir.ref<i32> {fir.bindc_name = "d"},
    ! CHECK-SAME:  %[[VAL_4:.*]]: !fir.ref<i32> {fir.bindc_name = "e", fir.optional}) {
    subroutine dynamic_optional_weird(a, b, c, d, e)
      integer :: a, b, d
      integer, optional :: c, e
      ! a3, a4, a6, a8 statically missing. a5, a9 dynamically optional.
      print *, max(a1=a, a2=b, a5=c, a7=d, a9 = e)
    ! CHECK:  %[[VAL_10:.*]] = fir.load %[[VAL_0]] : !fir.ref<i32>
    ! CHECK:  %[[VAL_11:.*]] = fir.load %[[VAL_1]] : !fir.ref<i32>
    ! CHECK:  %[[VAL_12:.*]] = fir.is_present %[[VAL_2]] : (!fir.ref<i32>) -> i1
    ! CHECK:  %[[VAL_13:.*]] = fir.load %[[VAL_3]] : !fir.ref<i32>
    ! CHECK:  %[[VAL_14:.*]] = fir.is_present %[[VAL_4]] : (!fir.ref<i32>) -> i1
    ! CHECK:  %[[VAL_15:.*]] = arith.cmpi sgt, %[[VAL_10]], %[[VAL_11]] : i32
    ! CHECK:  %[[VAL_16:.*]] = arith.select %[[VAL_15]], %[[VAL_10]], %[[VAL_11]] : i32
    ! CHECK:  %[[VAL_17:.*]] = fir.if %[[VAL_12]] -> (i32) {
    ! CHECK:    %[[VAL_18:.*]] = fir.load %[[VAL_2]] : !fir.ref<i32>
    ! CHECK:    %[[VAL_19:.*]] = arith.cmpi sgt, %[[VAL_16]], %[[VAL_18]] : i32
    ! CHECK:    %[[VAL_20:.*]] = arith.select %[[VAL_19]], %[[VAL_16]], %[[VAL_18]] : i32
    ! CHECK:    fir.result %[[VAL_20]] : i32
    ! CHECK:  } else {
    ! CHECK:    fir.result %[[VAL_16]] : i32
    ! CHECK:  }
    ! CHECK:  %[[VAL_21:.*]] = arith.cmpi sgt, %[[VAL_17]], %[[VAL_13]] : i32
    ! CHECK:  %[[VAL_23:.*]] = arith.select %[[VAL_21]], %[[VAL_17]], %[[VAL_13]] : i32
    ! CHECK:  %[[VAL_24:.*]] = fir.if %[[VAL_14]] -> (i32) {
    ! CHECK:    %[[VAL_25:.*]] = fir.load %[[VAL_4]] : !fir.ref<i32>
    ! CHECK:    %[[VAL_26:.*]] = arith.cmpi sgt, %[[VAL_23]], %[[VAL_25]] : i32
    ! CHECK:    %[[VAL_27:.*]] = arith.select %[[VAL_26]], %[[VAL_23]], %[[VAL_25]] : i32
    ! CHECK:    fir.result %[[VAL_27]] : i32
    ! CHECK:  } else {
    ! CHECK:    fir.result %[[VAL_23]] : i32
    ! CHECK:  }
    ! CHECK:  fir.call @_FortranAioOutputInteger32(%{{.*}}, %[[VAL_24]]) : (!fir.ref<i8>, i32) -> i1
    end subroutine
    end module

      use :: max_test
      integer :: a(4) = [1,12,23, 34]
      integer :: b(4) = [31,22,13, 4]
      integer :: c(4) = [21,32,3, 14]
      call dynamic_optional(a, b)
      call dynamic_optional(a, b, c)
      call dynamic_optional_array_expr_scalar_optional(a, b)
      call dynamic_optional_array_expr_scalar_optional(a, b, c(2))
      call dynamic_optional_scalar(a(2), b(2))
      call dynamic_optional_scalar(a(2), b(2), c(2))
    end

max_stars_repo_name: colemancda/llvm-project

max_issues_repo_name: colemancda/llvm-project

max_forks_repo_name: colemancda/llvm-project


text:
// --------------------------Synchronized-------------------------- //

//    Memory 1    //    Memory 2    //
// cw0            // cw1            //
// cw1            // cw2            //
// cw2            // cw3            //
// cw3            // cw4            //
// cw4            // cw5            //
// cw5            // Speed          //
// Speed          // Delta          //
// Delta          // Stabilization  //
// Stabilization  // FlushReceivers //
// FlushReceivers // cw-            //

// Scheduler Start - Start
:cw0=0 goto1

// Scheduler Stop - cw0
:cw0=-1 goto1

// Flush Receivers - FlushReceivers
:TargetFrequency=2 :TargetFrequency=1 goto1

// Stabilization - Stabilization
:O1=:O :F1=:F :U1=:U :O2=:P :F2=:G :U2=:V :O3=:Q :F3=:H :U3=:W goto1

// Trilateration Delta - Delta
O=241868 P=14625000000 Q=128996 R=145121 M=1000000 T=9425000000 h=0.5
N=M-:O3 E=M-:O1 W=M-:O2 N*=N E*=E W*=W :DX=X-A :DY=Y-B :DZ=Z-C
A=X B=Y C=Z X=(E-W+P)/O Y=(E-N+T-Q*X)/R Z=:ZSg*(E-X^2-Y^2)^h goto2

// Speed - Speed
:vf=X*:fi+Y*:fj+Z*:fk :vu=X*:ui+Y*:uj+Z*:uk :vr=X*:ri+Y*:rj+Z*:rk
d=0.6 X=:DX/d Y=:DY/d Z=:DZ/d goto1

// --------------------------Synchronized-------------------------- //

// --------------------------Asynchronized------------------------- //

// Signal Stability Indicator
t=2.5^2 A=:F1-:O1 B=:U1-:O1 C=:F2-:O2 D=:U2-:O2 E=:F3-:O3 F=:U3-:O3
:Nav=(A^2<t)*(B^2<t)*(C^2<t)*(D^2<t)*(E^2<t)*(F^2<t) goto1

// Coordinate Stability Indicator
A=abs(:D4) B=abs(:D) C=abs(:D5) D=abs(:D2) E=abs(:D6) F=abs(:D3)
t=2.5 :Cord=(A<t)*(B<t)*(C<t)*(D<t)*(E<t)*(F<t)

// Orientation Stability Indicator
:Ori=abs(:ui*:fi+:uj*:fj+:uk*:fk)<0.087 goto1

// Trilateration - U
O=241868 P=14625000000 Q=128996 R=145121 M=1000000 T=9425000000 h=0.5
N=M-:U3 S=M-:BS E=M-:U1 W=M-:U2 N*=N S*=S E*=E W*=W K=:ZSg
:UX=(E-W+P)/O :UY=(E-N+T-Q*:UX)/R :UZ=K*(E-:UX^2-:UY^2)^h
goto2

// Trilateration - F
O=241868 P=14625000000 Q=128996 R=145121 M=1000000 T=9425000000 h=0.5
N=M-:F3 S=M-:BS E=M-:F1 W=M-:F2 N*=N S*=S E*=E W*=W K=:ZSg
:FX=(E-W+P)/O :FY=(E-N+T-Q*:FX)/R :FZ=K*(E-:FX^2-:FY^2)^h
goto2

// Trilateration - O
O=241868 P=14625000000 Q=128996 R=145121 M=1000000 T=9425000000 h=0.5
N=M-:O3 S=M-:BS E=M-:O1 W=M-:O2 N*=N S*=S E*=E W*=W K=:ZSg
:X=(E-W+P)/O :Y=(E-N+T-Q*:X)/R :Z=K*(E-:X^2-:Y^2)^h
:D=:UX-:X :D2=:UY-:Y :D3=:UZ-:Z :D4=:FX-:X :D5=:FY-:Y :D6=:FZ-:Z goto2

// Trilateration - Reflection
G=(66359-:X)^2+(-37279-:Y)^2 R=(1000000-:BS) n=71111 :ZSg=1/(:ZSg==0)
A=abs(R-(G+(n-abs(:Z))^2)^.5) B=abs(R-(G+(n+abs(:Z))^2)^.5)
:ZSg = B>A - A>B goto1

// Orientation - U and F
h=.5
X=:D Y=:D2 Z=:D3 N=(X^2+Y^2+Z^2)^h A=:D4 B=:D5 C=:D6 M=(A^2+B^2+C^2)^h
:ui=X/N :uj=Y/N :uk=Z/N :fi=A/M :fj=B/M :fk=C/M goto2

// Orientation - R
:ri=:fj*:uk-:fk*:uj :rj=:fk*:ui-:fi*:uk :rk=:fi*:uj-:fj*:ui goto1

// Direction - Distance
:DWPx=:WPCX-:X :DWPy=:WPCY-:Y :DWPz=:WPCZ-:Z h=.5
:DWPD=(:DWPx^2+:DWPx^2+:DWPx^2)^h goto1

// Direction - Projection
x=:DWPx y=:DWPy z=:DWPz :WPCR=R :WPCF=F :WPCU=U
R=x*:ri+y*:rj+z*:rk F=x*:fi+y*:fj+z*:fk U=x*:ui+y*:uj+z*:uk goto1

// Direction - Yaw
h=.5 b=57.296 i=8/3 j=1/3
R=:WPCR F=:WPCF a=(F^2+R^2)^h a=F/a k=R>0-0>R
c=i*(2-(2+2*a)^h)^h-j*(2-2*a)^h :Yaw=k*c*b goto2

// Direction - Pitch
h=.5 b=57.296 i=0.281
R=:WPCR F=:WPCF U=:WPCU  t=U/(F^2+R^2)^h k=t^2>1 m=1-2*(t<0)
t^=1-2*k t*=b/(1+i*t^2) :Pit=m*k*90+(k<1-k)*t goto2

// --------------------------Asynchronized------------------------- //

// -------------------------Display Drivers------------------------ //

// Display - Location
h="Signa :1  : Trope :   : Coord :"
:1=h+"\nX:"+(:X/10*10)+"\nY:"+(:Y/10*10)+"\nZ:"+(:Z/10*10) goto1

// Display - Delta
h="Signa :2  : Trope :   : Delta :"
:2=h+"\nDX:"+(:DX/10*10)+"\nDY:"+(:DY/10*10)+"\nDZ:"+(:DZ/10*10) goto1

// Display - Speed
S="Signa :3  : Trope :   : Speed : "+"\nFB:"+:vF
:3=S+"\nRL:"+:vR+"\nUD:"+:vU+"\n    m/s" goto1

// Display - WayPoint
a=1000 h="Signa :4  : Trope :   : Dest. :"+"\nDX:"+(:DWPx/a*a)+"\nDY:"
:4=h+(:DWPy/a*a)+"\nDZ:"+(:DWPz/a*a)+"\nD:"+(:DWPD/a*a)+" m" goto1

// Display - CrossHair Horizontal
s=" " c="×"
i=(:Yaw/10)/1000*1000 if i>4 then i=4 end if i<-4 then i=-4 endgotoi+7
:horcom=c+s+s+s+s+s+s+s+s :navcentered=0 goto2
:horcom=s+c+s+s+s+s+s+s+s goto2
:horcom=s+s+c+s+s+s+s+s+s goto2
:horcom=s+s+s+c+s+s+s+s+s goto2
p=c if :verctr then p="※" end  :horcom=s+s+s+s+p+s+s+s+s goto2
:horcom=s+s+s+s+s+c+s+s+s goto2
:horcom=s+s+s+s+s+s+c+s+s goto2
:horcom=s+s+s+s+s+s+s+c+s goto2
:horcom=s+s+s+s+s+s+s+s+c goto2

// Display - CrossHair Vertical
t="  ¦   :0 " b="0:   ¦   :0 " l="            " m="-    +    - "
j=(:Pit/10)/1000*1000 if j>2 then j=2 end if j<-2 then j=-2 endgotoj+5
c=" "+:horcom+"  " :0=t+l+l+m+l+c+b :verctr=0 goto2
c=" "+:horcom+"  " :0=t+l+l+m+c+l+b :verctr=0 goto2
c="-"+:horcom+"- " :0=t+l+l+c+l+l+b :verctr=1 goto2
c=" "+:horcom+"  " :0=t+l+c+m+l+l+b :verctr=0 goto2
c=" "+:horcom+"  " :0=t+c+l+m+l+l+b :verctr=0 goto2

// -------------------------Display Drivers------------------------ //


text:
! ###############################################################
! #                                                             #
! #                       VLIDORT_2p8p3                         #
! #                                                             #
! #  Vectorized LInearized Discrete Ordinate Radiative Transfer #
! #  -          --         -        -        -         -        #
! #                                                             #
! ###############################################################

! ###############################################################
! #                                                             #
! #  Authors :     Robert. J. D. Spurr (1)                      #
! #                Matt Christi                                 #
! #                                                             #
! #  Address (1) : RT Solutions, inc.                           #
! #                9 Channing Street                            #
! #                Cambridge, MA 02138, USA                     #
! #                                                             #
! #  Tel:          (617) 492 1183                               #
! #  Email :       rtsolutions@verizon.net                      #
! #                                                             #
! #  This Version :   VLIDORT_2p8p3                             #
! #  Release Date :   31 March 2021                             #
! #                                                             #
! #  Previous VLIDORT Versions under Standard GPL 3.0:          #
! #  ------------------------------------------------           #
! #                                                             #
! #      2.7   F90, released        August 2014                 #
! #      2.8   F90, released        May    2017                 #
! #      2.8.1 F90, released        August 2019                 #
! #      2.8.2 F90, limited release May    2020                 #
! #                                                             #
! #  Features Summary of Recent VLIDORT Versions:               #
! #  -------------------------------------------                #
! #                                                             #
! #      NEW: TOTAL COLUMN JACOBIANS         (2.4)              #
! #      NEW: BPDF Land-surface KERNELS      (2.4R)             #
! #      NEW: Thermal Emission Treatment     (2.4RT)            #
! #      Consolidated BRDF treatment         (2.4RTC)           #
! #      f77/f90 Release                     (2.5)              #
! #      External SS / New I/O Structures    (2.6)              #
! #                                                             #
! #      SURFACE-LEAVING / BRDF-SCALING      (2.7)              #
! #      TAYLOR Series / OMP THREADSAFE      (2.7)              #
! #      New Water-Leaving Treatment         (2.8)              #
! #      LBBF & BRDF-Telescoping, enabled    (2.8)              #
! #      Several Performance Enhancements    (2.8)              #
! #      Water-leaving coupled code          (2.8.1)            #
! #      Planetary problem, media properties (2.8.1)            #
! #      Doublet geometry post-processing    (2.8.2)            #
! #      Reduction zeroing, dynamic memory   (2.8.2)            #
! #                                                             #
! #  Features Summary of This VLIDORT Version                   #
! #  ----------------------------------------                   #
! #                                                             #
! #   2.8.3, released 31 March 2021.                            #
! #     ==> Green's function RT solutions (Nstokes = 1 or 3)    #
! #     ==> Sphericity Corrections using MS source terms        #
! #     ==> BRDF upgrades, including new snow reflectance       #
! #     ==> SLEAVE Upgrades, extended water-leaving treatment   #
! #                                                             #
! ###############################################################

! ###################################################################
! #                                                                 #
! # This is Version 2.8.3 of the VLIDORT_2p8 software library.      #
! # This library comes with the Standard GNU General Public License,#
! # Version 3.0, 29 June 2007. Please read this license carefully.  #
! #                                                                 #
! #      VLIDORT Copyright (c) 2003-2021.                           #
! #          Robert Spurr, RT Solutions, Inc.                       #
! #          9 Channing Street, Cambridge, MA 02138, USA.           #
! #                                                                 #
! # This file is part of VLIDORT_2p8p3 ( Version 2.8.3 )            #
! #                                                                 #
! # VLIDORT_2p8p3 is free software: you can redistribute it         #
! # and/or modify it under the terms of the Standard GNU GPL        #
! # (General Public License) as published by the Free Software      #
! # Foundation, either version 3.0 of the License, or any           #
! # later version.                                                  #
! #                                                                 #
! # VLIDORT_2p8p3 is distributed in the hope that it will be        #
! # useful, but WITHOUT ANY WARRANTY; without even the implied      #
! # warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR         #
! # PURPOSE. See the Standard GNU General Public License (GPL)      #
! # for more details.                                               #
! #                                                                 #
! # You should have received a copy of the Standard GNU General     #
! # Public License (GPL) Version 3.0, along with the VLIDORT_2p8p3  #
! # code package. If not, see <http://www.gnu.org/licenses/>.       #
! #                                                                 #
! ###################################################################

program VSLEAVEplus_Tester

!  This is the Version 2.8 driver. Created 9/18/16 from 2.7 Driver
!     R. Spurr. RT Solutions Inc.

!  Upgrade for Version 2.8.1, August 2019
!  ---------------------------------------

!  Upgrade for Version 2.8.3, March 2021
!  -------------------------------------

!  Module files for VSLEAVE and VLIDORT. Strict usage, Version 2.8 upwards

   USE VSLEAVE_SUP_AUX_m, Only : VSLEAVE_READ_ERROR
   USE VSLEAVE_SUP_MOD_m
   USE VSLEAVE_LINSUP_MOD_m

   USE VLIDORT_PARS_m
   USE VLIDORT_IO_DEFS_m
   USE VLIDORT_LIN_IO_DEFS_m

   USE VLIDORT_VSLEAVE_SUP_ACCESSORIES_m
   USE VLIDORT_VSLEAVE_LINSUP_ACCESSORIES_m

   USE VLIDORT_AUX_m,    Only : VLIDORT_READ_ERROR, VLIDORT_WRITE_STATUS
   USE VLIDORT_INPUTS_m, Only : VLIDORT_INPUT_MASTER, VLIDORT_BRDF_Sup_Init, VLIDORT_SS_Sup_Init
   USE VLIDORT_MASTERS_m

   USE VLIDORT_L_INPUTS_m, Only : VLIDORT_L_INPUT_MASTER, VLIDORT_BRDF_LinSup_Init, VLIDORT_SS_LinSup_Init
   USE VLIDORT_LCS_MASTERS_m

   IMPLICIT NONE

!  VLIDORT file inputs status structure

   TYPE(VLIDORT_Input_Exception_Handling) :: VLIDORT_InputStatus

!  VLIDORT input structures

   TYPE(VLIDORT_Fixed_Inputs)             :: VLIDORT_FixIn
   TYPE(VLIDORT_Modified_Inputs)          :: VLIDORT_ModIn

!  VLIDORT supplements i/o structure

   TYPE(VLIDORT_Sup_InOut)                :: VLIDORT_Sup

!  VLIDORT output structure

   TYPE(VLIDORT_Outputs)                  :: VLIDORT_Out

!  VLIDORT linearized input structures

   TYPE(VLIDORT_Fixed_LinInputs)          :: VLIDORT_LinFixIn
   TYPE(VLIDORT_Modified_LinInputs)       :: VLIDORT_LinModIn

!  VLIDORT linearized supplements i/o structure

   TYPE(VLIDORT_LinSup_InOut)             :: VLIDORT_LinSup

!  VLIDORT linearized output structure

   TYPE(VLIDORT_LinOutputs)               :: VLIDORT_LinOut

!  VSLEAVE supplement file inputs status structure

   TYPE(VSLEAVE_Input_Exception_Handling) :: VSLEAVE_Sup_InputStatus

!  VSLEAVE supplement input structures

   TYPE(VSLEAVE_Sup_Inputs)               :: VSLEAVE_Sup_In

!  VSLEAVE supplement linearized input structure

   TYPE(VSLEAVE_LinSup_Inputs)            :: VSLEAVE_LinSup_In

!  VSLEAVE supplement output structure

   TYPE(VSLEAVE_Sup_Outputs)              :: VSLEAVE_Sup_Out
   TYPE(VSLEAVE_Output_Exception_Handling):: VSLEAVE_Sup_OutputStatus

!  VSLEAVE supplement linearized output structure

   TYPE(VSLEAVE_LinSup_Outputs)           :: VSLEAVE_LinSup_Out

!  VSLEAVE supplement / VLIDORT VSLEAVE-related inputs consistency check status

   TYPE(VLIDORT_Exception_Handling)       :: VLIDORT_VSLEAVECheck_Status

!  Local Variables
!  ===============

   integer, parameter :: maxlambdas = 900
   integer, parameter :: maxlevels  = maxlayers + 1

   integer :: nlambdas, data_nheights, nlayers, ngreek_moments_input, ngeometries
   integer :: i, l, n, na, nf, v, w

   logical :: openfileflag
   logical :: do_fluorescence, do_waterleaving
   logical :: do_observation_geometry, do_doublet_geometry

   double precision :: aerwt, raywt, raysca, aersca, aertau, tottau, totsca
   double precision :: q24, g24, qgr, ggr, ang, aerang, aerssa, aerasy
   double precision :: cumaer1, cumaer2, aerprv1, aerprv2, aers_550, aert_550
   double precision :: rho_1, rho_2, rho_a, diff, temp, psurf, ray2
   double precision :: aerphas(maxmoments_input),aer550(maxlayers)

   double precision :: lambdas(maxlambdas), lambda_start
   double precision :: rayleigh_xsec(maxlambdas),rayleigh_depol(maxlambdas)

   double precision :: height(maxlevels),tempr(maxlevels)
   double precision :: data_pressures(maxlevels),data_heights(maxlevels)
   double precision :: data_temperatures(maxlevels), earth_radius

   double precision :: height_grid(0:maxlayers),layer_temperatures(maxlayers)
   double precision :: layer_aircolumns(maxlayers)

   double precision :: deltau_vert_input(maxlayers)
   double precision :: omega_total_input(maxlayers)
   double precision :: greekmat_total_input(0:maxmoments_input,maxlayers,16)

   !double precision :: fmatrix_up ( maxlayers, max_geometries, 6 )
   !double precision :: fmatrix_dn ( maxlayers, max_geometries, 6 )

   double precision :: lambertian_albedo

   character(len=7) :: c1
   character(len=5) :: c2

!  Saved results

   double precision :: I_base(maxlambdas,max_geometries),&
                       I_pert(max_surfacewfs+max_sleavewfs,maxlambdas,max_geometries),&
                       Jac   (max_surfacewfs+max_sleavewfs,maxlambdas,max_geometries)

!  Loschmidt's number (particles/cm3).

   double precision, parameter :: RHO_STANDARD = 2.68675D+19

!  STP values.

   double precision, parameter :: PZERO  = 1013.25D0, TZERO  = 273.15D0, PTZERO = PZERO / TZERO
   double precision, parameter :: RHO_ZERO = RHO_STANDARD * TZERO / PZERO, CONSTANT = 1.0D+05 * RHO_ZERO

!  CO2 PPMV mixing ratio (for the Rayleigh stuff)

   double precision, parameter :: CO2_PPMV_MIXRATIO = 390.0d0

!  Help variables

   integer :: g, j1, j2, jw, k, m, nsl, o1, uta, q, nmi, n_totalsurface_wfs
   double precision :: eps, epsfac, epsfac1, epsfac2
   double precision :: par(2)

   character(len=120) :: trace
   character(len=200) :: Vsleave_DataPath

!  1/31/21. Version 2.8.3, introduce do_debug_input control

   LOGICAL :: do_debug_input

!  Two TESTS, one for Land (SIF) and the other for water (Water-leaving)

   LOGICAL          :: DO_TEST1 = .true.
   LOGICAL          :: DO_TEST2 = .true.

!  TEST 1 - Full Calculcation for FLUORESCENCE
!  ===========================================

!  Avoid if not flagged - go to 87

   if ( .not. DO_TEST1 ) go to 87

!  Initialize by-hand input flag
!  1/31/21. Version 2.8.3. Local flag for debug input

   do_debug_input = .false.

!  Progress

   write(*,'(1x,a/)')' DOING TEST 1 - Full Calculcation for FLUORESCENCE'

!  Problem setup for TEST 1
!  ------------------------

!  Set VSLEAVE data path

   Vsleave_DataPath = 'vlidort_v_test/data'

!  Wavelengths

   lambda_start = 640.0d0
   nlambdas = 181
   do w = 1, nlambdas
      lambdas(w) = lambda_start + 1.0d0*dble(w-1)
   enddo

!  Get the pressures, temperatures and heights from FILE.

   OPEN(1,FILE=Trim(Vsleave_DataPath)//'/temp_psurf.prf_35',STATUS='OLD')
   READ(1,*)DATA_NHEIGHTS,psurf
   DO I = 1, DATA_NHEIGHTS
      READ(1,*)HEIGHT(I), TEMPR(I)
   ENDDO
   CLOSE(1)

!  Set the pressures

   DATA_PRESSURES(DATA_NHEIGHTS)=psurf
   DO I = DATA_NHEIGHTS-1, 1 ,-1
      DATA_PRESSURES(I) = DATA_PRESSURES(&
       I+1)*exp(-9.81*28.9/8314.0*(height(DATA_NHEIGHTS-I+1)-&
       height(DATA_NHEIGHTS-I))*1.0E3/2.0*(1.0/tempr(DATA_NHEIGHTS-I+1)&
       +1.0/tempr(DATA_NHEIGHTS-I)))
   ENDDO

!  Set the data heights

   DO I = 1, DATA_NHEIGHTS
      DATA_HEIGHTS(I)=HEIGHT(DATA_NHEIGHTS-I+1)
      DATA_TEMPERATURES(I)=TEMPR(DATA_NHEIGHTS-I+1)
   ENDDO

!  Set the height grid, number of layers, Layer temperatures and Aircolumns

   HEIGHT_GRID = 0.0d0
   NLAYERS = DATA_NHEIGHTS - 1
   HEIGHT_GRID(0) = DATA_HEIGHTS(1)
   DO I = 1, NLAYERS
      RHO_1 = DATA_PRESSURES(I) / DATA_TEMPERATURES(I)
      RHO_2 = DATA_PRESSURES(I+1) / DATA_TEMPERATURES(I+1)
      TEMP  = 0.5D0*(DATA_TEMPERATURES(I)+DATA_TEMPERATURES(I+1))
      RHO_A = 0.5D0 * CONSTANT * ( RHO_1 + RHO_2 )
      DIFF  = DATA_HEIGHTS(I) - DATA_HEIGHTS(I+1)
      HEIGHT_GRID(I) = DATA_HEIGHTS(I+1)
      LAYER_TEMPERATURES(I) = TEMP
      LAYER_AIRCOLUMNS(I)   = DIFF * RHO_A
   ENDDO

!  Rayleigh function

   call Rayleigh_function &
    ( MAXLAMBDAS, CO2_PPMV_MIXRATIO, &
      NLAMBDAS,   LAMBDAS, &
      RAYLEIGH_XSEC, RAYLEIGH_DEPOL )

!  Loading areas

   aers_550 = 0.01d0
   q24 = 2.0d0 * aers_550 / (height_grid(0)-height_grid(24))
   g24 = q24 /  (height_grid(0)-height_grid(24))

   aert_550 = 0.15d0
   qgr = 2.0d0 * aert_550 / (height_grid(24)-height_grid(nlayers))
   ggr = qgr / (height_grid(24)-height_grid(nlayers))

!  Basic layering for aerosol

   aerprv1 = 0.0d0
   aerprv2 = 0.0d0
   do n = 1, nlayers
      if ( n.le.24 ) then
         cumaer1 = g24 * (height_grid(0)-height_grid(n))
         aer550(n) = 0.5*(cumaer1+aerprv1)*(height_grid(n-1)-height_grid(n))
         aerprv1 = cumaer1
      else
         cumaer2 = ggr * (height_grid(24)-height_grid(n))
         aer550(n) = 0.5*(cumaer2+aerprv2)*(height_grid(n-1)-height_grid(n))
         aerprv2 = cumaer2
      endif
   enddo

!  Other aerosol inputs

   aerAng = 1.01d0
   aerssa = 0.99d0
   aerasy = 0.75d0
   ngreek_moments_input = 50
   do l = 1, ngreek_moments_input
      aerphas(L) = dble(2*L+1) * aerasy ** dble(l)
   enddo

!  Define a proxy

   nmi = ngreek_moments_input

!  Initialize Test 1
!  -----------------

!  Initialize error file

   OPENFILEFLAG = .false.

!  Initialize perturbation variables

   eps = 1.0d-04
   epsfac = 1.0d0 + eps
   epsfac1 = 1.0001d0
   epsfac2 = 1.001d0

!  Initialize output

   I_base = 0.0d0
   I_pert = 0.0d0
   Jac    = 0.0d0

!  Start wavelength loop
!  ---------------------

   !do w = 1, nlambdas
   do w = 1, nlambdas, 10
   !do w = 1,1

      !write(*,*)
      !write(*,*) '******************************************************'
      write(*,*) 'w = ',w
      !write(*,*) 'lambdas(w) = ',lambdas(w)

      !write(*,*)
      !write(*,*) 'Doing Baseline Calc'
      !write(*,*) '-------------------'

!  @@@@@@@@@@@@@@@@@@@@@@@@@@@
!  BASELINE / ANALYTIC SECTION
!  @@@@@@@@@@@@@@@@@@@@@@@@@@@

!  Read VSLEAVE inputs
!  -------------------

      CALL VSLEAVE_LIN_INPUTMASTER ( &
        'vlidort_v_test/V2p8p3_VSLEAVE_ReadInput_fluor.cfg', & ! Input
         VSLEAVE_Sup_In,         & ! Outputs
         VSLEAVE_LinSup_In,      & ! Outputs
         VSLEAVE_Sup_InputStatus ) ! Outputs

      IF ( VSLEAVE_Sup_InputStatus%SL_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VSLEAVE_READ_ERROR ( 'V2p8p3_VSLEAVE_ReadInput_lin_call_Test1.log', VSLEAVE_Sup_InputStatus )

!  Define some inputs not handled by VSLEAVE_LIN_INPUTMASTER

      VSLEAVE_Sup_In%SL_VSLEAVE_DATAPATH = Trim(Vsleave_DataPath)

!  mick fix 4/1/2015 - moved setting of VSLEAVE wavelength from below

      VSLEAVE_Sup_In%SL_FL_WAVELENGTH = lambdas(w) !must be in [nm]

!  Define local variables

      DO_FLUORESCENCE         = VSLEAVE_Sup_In%SL_DO_FLUORESCENCE
      DO_WATERLEAVING         = .not. VSLEAVE_Sup_In%SL_DO_FLUORESCENCE
      DO_OBSERVATION_GEOMETRY = VSLEAVE_Sup_In%SL_DO_USER_OBSGEOMS
      DO_DOUBLET_GEOMETRY     = VSLEAVE_Sup_In%SL_DO_DOUBLET_GEOMETRY
      nsl                     = VSLEAVE_LinSup_In%SL_N_SLEAVE_WFS

!  Define VSLEAVE outputs
!  ----------------------

!  Linearized VSLEAVE call - the output will be used for the
!                            linearized VLIDORT calculation

      CALL VSLEAVE_LIN_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_LinSup_In,      & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_LinSup_Out,     & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_LIN_MAINMASTER failed' ; go to 678
      endif

!  Read VLIDORT Main inputs
!  ------------------------

!  VLIDORT linearized control read input, abort if failed

      CALL VLIDORT_L_INPUT_MASTER ( &
        'vlidort_v_test/V2p8p3_vsleaveplus_tester_fluor.cfg', & ! Input
        VLIDORT_FixIn,      & ! Outputs
        VLIDORT_ModIn,      & ! Outputs
        VLIDORT_LinFixIn,   & ! Outputs
        VLIDORT_LinModIn,   & ! Outputs
        VLIDORT_InputStatus ) ! Outputs

      IF ( VLIDORT_InputStatus%TS_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_READ_ERROR ( 'V2p8p3_VLIDORT_ReadInput_lin_call_Test1.log', VLIDORT_InputStatus )

!  1/31/21. Version 2.8.3, New variables must be set by hand

      VLIDORT_FixIn%Cont%TS_ASYMTX_TOLERANCE = 1.0d-20

!  1/31/21. Version 2.8.3, Set the DO_MSSTS flag to generate output for the MS sphericity corrections

      VLIDORT_FixIn%Bool%TS_DO_MSSTS = .false.

!  There are very strict conditions for this specialist option, as follows :==>
!     1. Either DO_UPWELLING or DO_DNWELLING must be set, Not Both !!!!
!     2. DO_FULLRAD_MODE must be set
!     3. DO_OBSERVATION_GEOMETRY must be set, with N_USER_OBSGEOMS = 3
!     4. DO_FOCORR and DO_FOCORR_OUTGOING must both be set
!     5a. Upwelling  : N_USER_LEVELS = 1, and USER_LEVELS(1) = 0.0            [ TOA output only ]
!     5b. Downwelling: N_USER_LEVELS = 1, and USER_LEVELS(1) = Real(nlayers)  [ BOA output only ]
!  These checks have been implemented inside VLIDORT.

!  1/31/21. Version 2.8.3, Set the NSTOKES2 FOURIER0 flag by hand

      VLIDORT_FixIn%Bool%TS_DO_FOURIER0_NSTOKES2 = .true.

!  Redefine some VLIDORT input variables (override VLIDORT config file)
!  mick fix 4/1/2015 - added this wavelength diagnostic (just a dummy here)

      VLIDORT_FixIn%Optical%TS_ATMOS_WAVELENGTH = lambdas(w)/1.0D3 !must be in [um]

!  Define some input variables

      VLIDORT_LinFixIn%Cont%TS_N_SURFACE_WFS = 1
      VLIDORT_LinFixIn%Cont%TS_N_SLEAVE_WFS  = VSLEAVE_LinSup_In%SL_N_SLEAVE_WFS

!  Initialize some input variables not handled by VLIDORT_L_INPUT_MASTER

      !CALL VLIDORT_BRDF_Sup_Init ( VLIDORT_Sup )
      !CALL VLIDORT_SS_Sup_Init ( VLIDORT_Sup )
      !CALL VLIDORT_BRDF_LinSup_Init ( VLIDORT_LinSup )
      !CALL VLIDORT_SS_LinSup_Init ( VLIDORT_LinSup )

!  IMPORTANT - Check compatibility of VSLEAVE and VLIDORT Main inputs
!  ------------------------------------------------------------------

      CALL VLIDORT_VSLEAVE_INPUT_CHECK ( &
        VSLEAVE_Sup_In,             & ! Inputs
        VLIDORT_FixIn,              & ! Inputs
        VLIDORT_ModIn,              & ! Inputs
        VLIDORT_VSLEAVECheck_Status ) ! Outputs

!  Exception handling

      IF ( VLIDORT_VSLEAVECheck_Status%TS_STATUS_INPUTCHECK .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_VSLEAVE_INPUT_CHECK_ERROR ( &
           'V2p8p3_VLIDORT_VSLEAVEcheck_lin_call_Test1.log', VLIDORT_VSLEAVECheck_Status )

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std & lin)

      VLIDORT_LinModIn%MCont%TS_DO_SLEAVE_WFS = .TRUE.
      CALL SET_VLIDORT_L_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out,  VSLEAVE_LinSup_Out, & !Inputs
        VLIDORT_FixIn,    VLIDORT_ModIn,      & !Inputs
        VLIDORT_LinFixIn, VLIDORT_LinModIn,   & !Inputs
        VLIDORT_Sup, VLIDORT_LinSup )           !Outputs

      N_TOTALSURFACE_WFS = VLIDORT_LinFixIn%Cont%TS_N_SURFACE_WFS + &
                           VLIDORT_LinFixIn%Cont%TS_N_SLEAVE_WFS

!  Optical properties
!  ------------------

!  Angstrom and ray2mom

      ray2 = ( 1.0d0 - rayleigh_depol(w) ) / ( 2.0d0 + rayleigh_depol(w) )
      ang  = ( lambdas(w) / 550.0d0 ) ** (-aerAng)

!      write(*,*)w,lambdas(w),ang,rayleigh_depol(w)

!  Make VLIDORT optical properties

      do n = 1, nlayers
         aertau = aer550(n) * ang
         aersca = aertau * aerssa
         raysca = LAYER_AIRCOLUMNS(n) * rayleigh_xsec(w)
         tottau = raysca + aertau
         totsca = raysca + aersca
         deltau_vert_input(n) = tottau
         omega_total_input(n) = totsca / tottau
         raywt = raysca / totsca ; aerwt = 1.0d0 - raywt
         greekmat_total_input(0,n,1) = 1.0d0
         greekmat_total_input(1,n,1) = aerwt * aerphas(1)
         greekmat_total_input(2,n,1) = aerwt * aerphas(2) + raywt * ray2
         do L = 3, ngreek_moments_input
            greekmat_total_input(L,n,1) = aerwt * aerphas(L)
         enddo
      enddo

      lambertian_albedo = 0.02d0

!  Save baseline parameters

      par(1) = lambertian_albedo
      par(2) = VSLEAVE_Sup_In%SL_FL_Amplitude755

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo    = lambertian_albedo

!  Linearized VLIDORT call
!  ----------------------
!  1/31/21. Version 2.8.3. Need to set the debug_input flag

      do_debug_input = .false.
      CALL VLIDORT_LCS_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup,   &
        VLIDORT_Out,   &
        VLIDORT_LinFixIn, &
        VLIDORT_LinModIn, &
        VLIDORT_LinSup,   &
        VLIDORT_LinOut )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save baseline results
!  ---------------------

      uta = 1; o1 = 1
      do v = 1, VSLEAVE_Sup_In%SL_NBEAMS
        I_base(w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
        do q = 1,n_totalsurface_wfs
          Jac(q,w,v) = VLIDORT_LinOut%Surf%TS_surfacewf(q,uta,v,o1,upidx)
        enddo
      enddo

!  @@@@@@@@@@@@@@@@@@@@@@@@@
!  FINITE DIFFERENCE SECTION
!  @@@@@@@@@@@@@@@@@@@@@@@@@

      !write(*,*)
      !write(*,*) 'Doing FD Calc'
      !write(*,*) '-------------'

!  Read the Vsleave and VLIDORT configuration files (again)
!  --------------------------------------------------------

!  Get the VSLEAVE inputs

      CALL VSLEAVE_INPUTMASTER ( &
        'vlidort_v_test/V2p8p3_VSLEAVE_ReadInput_fluor.cfg', & ! Input
        VSLEAVE_Sup_In,                         & ! Outputs
        VSLEAVE_Sup_InputStatus )                 ! Outputs

      IF ( VSLEAVE_Sup_InputStatus%SL_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VSLEAVE_READ_ERROR ( 'V2p8p3_VSLEAVE_ReadInput_fd_call_Test1.log', VSLEAVE_Sup_InputStatus )

!  Define some inputs not handled by VSLEAVE_INPUTMASTER

      VSLEAVE_Sup_In%SL_VSLEAVE_DATAPATH = Trim(Vsleave_DataPath)

!  mick fix 4/1/2015 - set VSLEAVE wavelength here

      VSLEAVE_Sup_In%SL_FL_WAVELENGTH = lambdas(w) !must be in [nm]

!  VLIDORT control Read input, abort if failed

      CALL VLIDORT_INPUT_MASTER ( &
        'vlidort_v_test/V2p8p3_vsleaveplus_tester_fluor.cfg', & ! Input
        VLIDORT_FixIn,      & ! Outputs
        VLIDORT_ModIn,      & ! Outputs
        VLIDORT_InputStatus ) ! Outputs

      IF ( VLIDORT_InputStatus%TS_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_READ_ERROR ( 'V2p8p3_VLIDORT_ReadInput_fd_call_Test1.log', VLIDORT_InputStatus )

!  1/31/21. Version 2.8.3, New variables must be set by hand

      VLIDORT_FixIn%Cont%TS_ASYMTX_TOLERANCE = 1.0d-20

!  1/31/21. Version 2.8.3, Set the DO_MSSTS flag to generate output for the MS sphericity corrections

      VLIDORT_FixIn%Bool%TS_DO_MSSTS = .false.

!  There are very strict conditions for this specialist option, as follows :==>
!     1. Either DO_UPWELLING or DO_DNWELLING must be set, Not Both !!!!
!     2. DO_FULLRAD_MODE must be set
!     3. DO_OBSERVATION_GEOMETRY must be set, with N_USER_OBSGEOMS = 3
!     4. DO_FOCORR and DO_FOCORR_OUTGOING must both be set
!     5a. Upwelling  : N_USER_LEVELS = 1, and USER_LEVELS(1) = 0.0            [ TOA output only ]
!     5b. Downwelling: N_USER_LEVELS = 1, and USER_LEVELS(1) = Real(nlayers)  [ BOA output only ]
!  These checks have been implemented inside VLIDORT.

!  1/31/21. Version 2.8.3, Set the NSTOKES2 FOURIER0 flag by hand

      VLIDORT_FixIn%Bool%TS_DO_FOURIER0_NSTOKES2 = .true.

!  mick fix 4/1/2015 - added wavelength diagnostic (just a dummy here)
!  Redefine some VLIDORT input variables (override VLIDORT config file)

      VLIDORT_FixIn%Optical%TS_ATMOS_WAVELENGTH = lambdas(w)/1.0D3 !must be in [um]

!  Initialize some input variables not handled by VLIDORT_INPUT_MASTER

      !CALL VLIDORT_BRDF_Sup_Init ( VLIDORT_Sup )
      !CALL VLIDORT_SS_Sup_Init ( VLIDORT_Sup )

!   IMPORTANT - Check compatibility of VSLEAVE and Main VLIDORT inputs
!   ------------------------------------------------------------------

      CALL VLIDORT_VSLEAVE_INPUT_CHECK ( &
        VSLEAVE_Sup_In,             & ! Inputs
        VLIDORT_FixIn,              & ! Inputs
        VLIDORT_ModIn,              & ! Inputs
        VLIDORT_VSLEAVECheck_Status ) ! Outputs

      IF ( VLIDORT_VSLEAVECheck_Status%TS_STATUS_INPUTCHECK .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_VSLEAVE_INPUT_CHECK_ERROR ( &
           'V2p8p3_VLIDORT_VSLEAVEcheck_fd_call_Test1.log', VLIDORT_VSLEAVECheck_Status )

!  Finite diff jacobian #1: wrt FL_Amplitude755
!  --------------------------------------------

!  Perturb FL_Amplitude755

      VSLEAVE_Sup_In%SL_FL_Amplitude755 = par(2)*epsfac

!  VSLEAVE call - the output will be used for the VLIDORT calculation

      CALL VSLEAVE_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_MAINMASTER failed, perturbation # 1' ; go to 678
      endif

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std only)

      CALL SET_VLIDORT_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out, VLIDORT_FixIn, VLIDORT_ModIn, & !Inputs
        VLIDORT_Sup )                                    !Outputs

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo                        = lambertian_albedo

!  VLIDORT call
!  1/31/21. Version 2.8.3. Need to set the debug_input flag

      do_debug_input = .false.
      CALL VLIDORT_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup, &
        VLIDORT_Out )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save perturbed results

      q = 2 !for FL_Amplitude755 (follows Lambertian albedo jacobian in array)
      do v = 1, VSLEAVE_Sup_In%SL_NBEAMS
        I_pert(q,w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
      enddo

!      write(*,*)
!      write(*,*) 'q = ',q
!      write(*,51)w,lambdas(w),' I_pert(q,w,v)                   = ',&
!                 (I_pert(q,w,v) ,v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' Jac(q,w,v)                      = ',&
!                 (Jac(q,w,v)    ,v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' par(q)*Jac(q,w,v)               = ',&
!                 (par(q)*Jac(q,w,v),v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' (I_pert(q,w,v)-I_base(w,v))/eps = ',&
!                 ((I_pert(q,w,v)-I_base(w,v))/eps,v=1,VSLEAVE_Sup_In%SL_NBEAMS)

!  Reset FL_Amplitude755

      VSLEAVE_Sup_In%SL_FL_Amplitude755 = par(2)

!  Finite diff jacobian #2: wrt Lambertian albedo
!  ----------------------------------------------

!  VSLEAVE call - the output will be used for the VLIDORT calculation

      CALL VSLEAVE_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_MAINMASTER failed, perturbation # 2' ; go to 678
      endif

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std only)

      CALL SET_VLIDORT_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out, VLIDORT_FixIn, VLIDORT_ModIn, & !Inputs
        VLIDORT_Sup )                                    !Outputs

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo    = lambertian_albedo

!  Perturb Lambertian albedo

      VLIDORT_FixIn%Optical%TS_LAMBERTIAN_ALBEDO  = par(1)*epsfac

!  VLIDORT call
!  1/31/21. Version 2.8.3. Need to set the debug_input flag

      do_debug_input = .false.
      CALL VLIDORT_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup, &
        VLIDORT_Out )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save perturbed results

      q = 1 !for lambertian albedo jacobian (first in array)
      do v = 1, VSLEAVE_Sup_In%SL_NBEAMS
        I_pert(q,w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
      enddo

!      write(*,*)
!      write(*,*) 'q = ',q
!      write(*,51)w,lambdas(w),' I_pert(q,w,v)                   = ',&
!                 (I_pert(q,w,v) ,v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' Jac(q,w,v)                      = ',&
!                 (Jac(q,w,v)    ,v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' par(q)*Jac(q,w,v)               = ',&
!                 (par(q)*Jac(q,w,v),v=1,VSLEAVE_Sup_In%SL_NBEAMS)
!      write(*,51)w,lambdas(w),' (I_pert(q,w,v)-I_base(w,v))/eps = ',&
!                 ((I_pert(q,w,v)-I_base(w,v))/eps,v=1,VSLEAVE_Sup_In%SL_NBEAMS)

51 format(i4,f10.4,a,70e17.10)

!  End wavelength loop

   enddo

!  Write results, Test 1
!  ---------------------

!  Headers

   if (do_observation_geometry) then
     c1 = 'obsgeo'
   else if (do_doublet_geometry) then
     c1 = 'doublet'
   else
     c1 = 'lattice'
   endif

   if (do_fluorescence) then
     c2 = 'fluor'
   else
     c2 = 'clean'
   end if

!  Open file, write results

   open(unit=777,file='vlidort_v_test/results_vsleaveplus_tester_' &
        // trim(c1) // '_' // c2 // '.all',status='replace',action='write')

   !do w = 1, nlambdas
   do w = 1, nlambdas, 10
   !do w = 1,1
      write(777,52) w,lambdas(w),&
        ( I_base(w,v),( par(q)*Jac(q,w,v),(I_pert(q,w,v)-I_base(w,v))/eps,&
          q=1,n_totalsurface_wfs ),v=1,VSLEAVE_Sup_In%SL_NBEAMS )
   enddo
52 format(i4,f10.4,70e12.4)

   close(777)

!  Finish Test 1

   write(*,*)
   write(*,*) ' TEST 1 finished successfully'

! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!  continuation point for Test 2

87 continue

!  Avoid if not flagged - go to 89

   if ( .not. DO_TEST2 ) then
     go to 89
   else
     write(*,*)
   endif

!  TEST 2 - Full Calculcation for WATER-LEAVING
!  ============================================

!  Initialize by-hand input flag
!  1/31/21. Version 2.8.3. Local flag for debug input

   do_debug_input = .false.

!  5/24/21. Select geometry

   do_doublet_geometry     = .false.
   do_observation_geometry = .false.

!  Headers

   if ( do_observation_geometry ) then
     c1 = 'obsgeom'
   else if ( do_doublet_geometry ) then
     c1 = 'doublet'
   else
     c1 = 'lattice'
   endif

!  Progress

   write(*,'(1x,a/)')' DOING TEST 2 - Full Calculcation for WATER-LEAVING'

!  Problem setup for TEST 2
!  ------------------------

!  Set VSLEAVE data path

   Vsleave_DataPath = 'vlidort_v_test/data'

!  Wavelengths

   lambda_start = 400.0d0
   nlambdas     = 21
   do w = 1, nlambdas
      lambdas(w) = lambda_start + 1.0d0*dble(w-1)
   enddo

!  Get the pressures, temperatures and heights from FILE.

   OPEN(1,FILE=Trim(Vsleave_DataPath)//'/temp_psurf.prf_35',STATUS='OLD')
   READ(1,*)DATA_NHEIGHTS,psurf
   DO I = 1, DATA_NHEIGHTS
       READ(1,*)HEIGHT(I), TEMPR(I)
   ENDDO
   CLOSE(1)

!  Set the pressures

   DATA_PRESSURES(DATA_NHEIGHTS)=psurf
   DO I = DATA_NHEIGHTS-1, 1 ,-1
      DATA_PRESSURES(I) = DATA_PRESSURES(&
       I+1)*exp(-9.81*28.9/8314.0*(height(DATA_NHEIGHTS-I+1)-&
       height(DATA_NHEIGHTS-I))*1.0E3/2.0*(1.0/tempr(DATA_NHEIGHTS-I+1)&
       +1.0/tempr(DATA_NHEIGHTS-I)))
   ENDDO

!  Set the data heights

   DO I = 1, DATA_NHEIGHTS
      DATA_HEIGHTS(I)=HEIGHT(DATA_NHEIGHTS-I+1)
      DATA_TEMPERATURES(I)=TEMPR(DATA_NHEIGHTS-I+1)
   ENDDO

!  Set the height grid, number of layers, Layer temperatures and Aircolumns

   HEIGHT_GRID = 0.0d0
   NLAYERS = DATA_NHEIGHTS - 1
   HEIGHT_GRID(0) = DATA_HEIGHTS(1)
   DO I = 1, NLAYERS
      RHO_1 = DATA_PRESSURES(I) / DATA_TEMPERATURES(I)
      RHO_2 = DATA_PRESSURES(I+1) / DATA_TEMPERATURES(I+1)
      TEMP  = 0.5D0*(DATA_TEMPERATURES(I)+DATA_TEMPERATURES(I+1))
      RHO_A = 0.5D0 * CONSTANT * ( RHO_1 + RHO_2 )
      DIFF  = DATA_HEIGHTS(I) - DATA_HEIGHTS(I+1)
      HEIGHT_GRID(I) = DATA_HEIGHTS(I+1)
      LAYER_TEMPERATURES(I) = TEMP
      LAYER_AIRCOLUMNS(I)   = DIFF * RHO_A
   ENDDO

!  Rayleigh function

   call Rayleigh_function &
    ( MAXLAMBDAS, CO2_PPMV_MIXRATIO, &
      NLAMBDAS,   LAMBDAS, &
      RAYLEIGH_XSEC, RAYLEIGH_DEPOL )

!  Loading areas

   aers_550 = 0.01d0
   q24 = 2.0d0 * aers_550 / (height_grid(0)-height_grid(24))
   g24 = q24 /  (height_grid(0)-height_grid(24))

   aert_550 = 0.15d0
   qgr = 2.0d0 * aert_550 / (height_grid(24)-height_grid(nlayers))
   ggr = qgr / (height_grid(24)-height_grid(nlayers))

!  Basic layering for aerosol

   aerprv1 = 0.0d0
   aerprv2 = 0.0d0
   do n = 1, nlayers
      if ( n.le.24 ) then
         cumaer1 = g24 * (height_grid(0)-height_grid(n))
         aer550(n) = 0.5*(cumaer1+aerprv1)*(height_grid(n-1)-height_grid(n))
         aerprv1 = cumaer1
      else
         cumaer2 = ggr * (height_grid(24)-height_grid(n))
         aer550(n) = 0.5*(cumaer2+aerprv2)*(height_grid(n-1)-height_grid(n))
         aerprv2 = cumaer2
      endif
   enddo

!  Other aerosol inputs

   aerAng = 1.01d0
   aerssa = 0.99d0
   aerasy = 0.75d0
   ngreek_moments_input = 50
   do l = 1, ngreek_moments_input
      aerphas(L) = dble(2*L+1) * aerasy ** dble(l)
   enddo

!  Define a proxy

   nmi = ngreek_moments_input

!  Initialize Test 2
!  -----------------

!  Initialize error file

   OPENFILEFLAG = .false.

!  Initialize perturbation variables

   eps = 1.0d-03
   epsfac = 1.0d0 + eps
   epsfac1 = epsfac
   epsfac2 = epsfac

!  Initialize output

   I_base = 0.0d0
   I_pert = 0.0d0
   Jac    = 0.0d0

!  Start wavelength loop
!  ---------------------

   do w = 1, nlambdas
   !do w = 1, 1

!  Progress

      write(*,*) 'w = ',w
      !write(*,*) 'lambdas(w) = ',lambdas(w)
      !write(*,'(A,I5)') '  Test 2: Doing wavelength # w = ',w

!  @@@@@@@@@@@@@@@@@@@@@@@@@@@
!  BASELINE / ANALYTIC SECTION
!  @@@@@@@@@@@@@@@@@@@@@@@@@@@

!  Read VSLEAVE inputs
!  -------------------

      CALL VSLEAVE_LIN_INPUTMASTER ( &
        'vlidort_v_test/V2p8p3_VSLEAVE_ReadInput_water.cfg', & ! Input
         VSLEAVE_Sup_In,         & ! Outputs
         VSLEAVE_LinSup_In,      & ! Outputs
         VSLEAVE_Sup_InputStatus ) ! Outputs

      IF ( VSLEAVE_Sup_InputStatus%SL_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VSLEAVE_READ_ERROR ( 'V2p8p3_VSLEAVE_ReadInput_lin_call_Test2.log', VSLEAVE_Sup_InputStatus )

!  Define some inputs not handled by VSLEAVE_LIN_INPUTMASTER

      VSLEAVE_Sup_In%SL_VSLEAVE_DATAPATH = Trim(Vsleave_DataPath)

!  Redefine VSLEAVE input wavelength (override VSLEAVE config file)

      VSLEAVE_Sup_In%SL_WAVELENGTH = lambdas(w)/1.0D3 !in Microns

!  Define local variables

      DO_FLUORESCENCE         = VSLEAVE_Sup_In%SL_DO_FLUORESCENCE
      DO_WATERLEAVING         = .not. VSLEAVE_Sup_In%SL_DO_FLUORESCENCE
      DO_OBSERVATION_GEOMETRY = VSLEAVE_Sup_In%SL_DO_USER_OBSGEOMS
      DO_DOUBLET_GEOMETRY     = VSLEAVE_Sup_In%SL_DO_DOUBLET_GEOMETRY
      nsl                     = VSLEAVE_LinSup_In%SL_N_SLEAVE_WFS

      if ( do_doublet_geometry ) then
         ngeometries = VSLEAVE_Sup_In%SL_NBEAMS * VSLEAVE_Sup_In%SL_N_USER_STREAMS
      else if ( do_observation_geometry ) then
         ngeometries = VSLEAVE_Sup_In%SL_NBEAMS
      else
         ngeometries = VSLEAVE_Sup_In%SL_NBEAMS * VSLEAVE_Sup_In%SL_N_USER_RELAZMS * VSLEAVE_Sup_In%SL_N_USER_STREAMS
      endif

!  Calculate VSLEAVE outputs
!  -------------------------

!  Linearized VSLEAVE call - the output will be used for the
!                            linearized VLIDORT calculation

      CALL VSLEAVE_LIN_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_LinSup_In,      & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_LinSup_Out,     & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_LIN_MAINMASTER failed' ; go to 678
      endif

!  Read VLIDORT Main inputs
!  ------------------------

!  VLIDORT linearized control read input, abort if failed

      CALL VLIDORT_L_INPUT_MASTER ( &
        'vlidort_v_test/V2p8p3_vsleaveplus_tester_water.cfg', & ! Input
        VLIDORT_FixIn,      & ! Outputs
        VLIDORT_ModIn,      & ! Outputs
        VLIDORT_LinFixIn,   & ! Outputs
        VLIDORT_LinModIn,   & ! Outputs
        VLIDORT_InputStatus ) ! Outputs

      IF ( VLIDORT_InputStatus%TS_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_READ_ERROR ( 'V2p8p3_VLIDORT_ReadInput_lin_call_Test2.log', VLIDORT_InputStatus )

!  1/31/21. Version 2.8.3, New variables must be set by hand

      VLIDORT_FixIn%Cont%TS_ASYMTX_TOLERANCE = 1.0d-20

!  1/31/21. Version 2.8.3, Set the DO_MSSTS flag to generate output for the MS sphericity corrections

      VLIDORT_FixIn%Bool%TS_DO_MSSTS = .false.

!  There are very strict conditions for this specialist option, as follows :==>
!     1. Either DO_UPWELLING or DO_DNWELLING must be set, Not Both !!!!
!     2. DO_FULLRAD_MODE must be set
!     3. DO_OBSERVATION_GEOMETRY must be set, with N_USER_OBSGEOMS = 3
!     4. DO_FOCORR and DO_FOCORR_OUTGOING must both be set
!     5a. Upwelling  : N_USER_LEVELS = 1, and USER_LEVELS(1) = 0.0            [ TOA output only ]
!     5b. Downwelling: N_USER_LEVELS = 1, and USER_LEVELS(1) = Real(nlayers)  [ BOA output only ]
!  These checks have been implemented inside VLIDORT.

!  1/31/21. Version 2.8.3, Set the NSTOKES2 FOURIER0 flag by hand

      VLIDORT_FixIn%Bool%TS_DO_FOURIER0_NSTOKES2 = .true.

!  Redefine some VLIDORT input variables (override VLIDORT config file)

      VLIDORT_FixIn%Optical%TS_ATMOS_WAVELENGTH = lambdas(w)/1.0D3 !must be in [um]

!  Define some inputs not handled by LIDORT_L_INPUT_MASTER

      VLIDORT_LinFixIn%Cont%TS_N_SURFACE_WFS = 0
      VLIDORT_LinFixIn%Cont%TS_N_SLEAVE_WFS  = VSLEAVE_LinSup_In%SL_N_SLEAVE_WFS

!  IMPORTANT - Check compatibility of VSLEAVE and VLIDORT Main inputs
!  ------------------------------------------------------------------

      CALL VLIDORT_VSLEAVE_INPUT_CHECK ( &
        VSLEAVE_Sup_In,             & ! Inputs
        VLIDORT_FixIn,              & ! Inputs
        VLIDORT_ModIn,              & ! Inputs
        VLIDORT_VSLEAVECheck_Status ) ! Outputs

!  Exception handling

      IF ( VLIDORT_VSLEAVECheck_Status%TS_STATUS_INPUTCHECK .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_VSLEAVE_INPUT_CHECK_ERROR ( &
           'V2p8p3_VLIDORT_VSLEAVEcheck_lin_call_Test2.log', VLIDORT_VSLEAVECheck_Status )

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std & lin)

      VLIDORT_LinModIn%MCont%TS_DO_SLEAVE_WFS = .TRUE.
      CALL SET_VLIDORT_L_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out,  VSLEAVE_LinSup_Out, & !Inputs
        VLIDORT_FixIn,    VLIDORT_ModIn,      & !Inputs
        VLIDORT_LinFixIn, VLIDORT_LinModIn,   & !Inputs
        VLIDORT_Sup, VLIDORT_LinSup )           !Outputs

      N_TOTALSURFACE_WFS = VLIDORT_LinFixIn%Cont%TS_N_SURFACE_WFS + &
                           VLIDORT_LinFixIn%Cont%TS_N_SLEAVE_WFS

!  Optical properties
!  ------------------

!  Angstrom and ray2mom

      ray2 = ( 1.0d0 - rayleigh_depol(w) ) / ( 2.0d0 + rayleigh_depol(w) )
      ang  = ( lambdas(w) / 550.0d0 ) ** (-aerAng)

      !write(*,*)w,lambdas(w),ang,rayleigh_depol(w)

!  Make VLIDORT optical properties

      do n = 1, nlayers
         aertau = aer550(n) * ang
         aersca = aertau * aerssa
         raysca = LAYER_AIRCOLUMNS(n) * rayleigh_xsec(w)
         tottau = raysca + aertau
         totsca = raysca + aersca
         deltau_vert_input(n) = tottau
         omega_total_input(n) = totsca / tottau
         raywt = raysca / totsca ; aerwt = 1.0d0 - raywt
         greekmat_total_input(0,n,1) = 1.0d0
         greekmat_total_input(1,n,1) = aerwt * aerphas(1)
         greekmat_total_input(2,n,1) = aerwt * aerphas(2) + raywt * ray2
         do L = 3, ngreek_moments_input
            greekmat_total_input(L,n,1) = aerwt * aerphas(L)
         enddo
      enddo

!stop 'oprop stop'

!  No albedo here

      lambertian_albedo = 0.0d0

!  Save baseline parameters

      par(1) = VSLEAVE_Sup_In%SL_CHLORCONC
      par(2) = VSLEAVE_Sup_In%SL_WINDSPEED

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo                        = lambertian_albedo

!  Linearized VLIDORT call
!  ----------------------
!  1/31/21. Version 2.8.3. Need to set the debug_input flag
!  1/31/21. Version 2.8.3. LCS here, works equally with LPS

      do_debug_input = .false.
      CALL VLIDORT_LCS_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup,   &
        VLIDORT_Out,   &
        VLIDORT_LinFixIn, &
        VLIDORT_LinModIn, &
        VLIDORT_LinSup,   &
        VLIDORT_LinOut )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save baseline results
!  ---------------------

      uta = 1; o1 = 1
      do v = 1, ngeometries
        I_base(w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
        do q = 1,n_totalsurface_wfs
          Jac(q,w,v) = VLIDORT_LinOut%Surf%TS_surfacewf(q,uta,v,o1,upidx)
        enddo
      enddo

!  @@@@@@@@@@@@@@@@@@@@@@@@@
!  FINITE DIFFERENCE SECTION
!  @@@@@@@@@@@@@@@@@@@@@@@@@

      !write(*,*)
      !write(*,*) 'Doing FD Calc'
      !write(*,*) '-------------'

!  Read the Vsleave and VLIDORT configuration files (again)
!  --------------------------------------------------------

!  Get the VSLEAVE inputs

      CALL VSLEAVE_INPUTMASTER ( &
        'vlidort_v_test/V2p8p3_VSLEAVE_ReadInput_water.cfg', & ! Input
        VSLEAVE_Sup_In,                         & ! Outputs
        VSLEAVE_Sup_InputStatus )                 ! Outputs

      IF ( VSLEAVE_Sup_InputStatus%SL_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VSLEAVE_READ_ERROR ( 'V2p8p3_VSLEAVE_ReadInput_fd_call_Test2.log', VSLEAVE_Sup_InputStatus )

!  Define some inputs not handled by VSLEAVE_INPUTMASTER

      VSLEAVE_Sup_In%SL_VSLEAVE_DATAPATH = Trim(Vsleave_DataPath)

!  Redefine VSLEAVE input wavelength (override VSLEAVE config file)

      VSLEAVE_Sup_In%SL_WAVELENGTH = lambdas(w)/1.0D3 !in Microns

!  VLIDORT control Read input, abort if failed

      CALL VLIDORT_INPUT_MASTER ( &
        'vlidort_v_test/V2p8p3_vsleaveplus_tester_water.cfg', & ! Input
        VLIDORT_FixIn,      & ! Outputs
        VLIDORT_ModIn,      & ! Outputs
        VLIDORT_InputStatus ) ! Outputs

      IF ( VLIDORT_InputStatus%TS_STATUS_INPUTREAD .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_READ_ERROR ( 'V2p8p3_VLIDORT_ReadInput_fd_call_Test2.log', VLIDORT_InputStatus )

!  1/31/21. Version 2.8.3, New variables must be set by hand

      VLIDORT_FixIn%Cont%TS_ASYMTX_TOLERANCE = 1.0d-20

!  1/31/21. Version 2.8.3, Set the DO_MSSTS flag to generate output for the MS sphericity corrections

      VLIDORT_FixIn%Bool%TS_DO_MSSTS = .false.

!  There are very strict conditions for this specialist option, as follows :==>
!     1. Either DO_UPWELLING or DO_DNWELLING must be set, Not Both !!!!
!     2. DO_FULLRAD_MODE must be set
!     3. DO_OBSERVATION_GEOMETRY must be set, with N_USER_OBSGEOMS = 3
!     4. DO_FOCORR and DO_FOCORR_OUTGOING must both be set
!     5a. Upwelling  : N_USER_LEVELS = 1, and USER_LEVELS(1) = 0.0            [ TOA output only ]
!     5b. Downwelling: N_USER_LEVELS = 1, and USER_LEVELS(1) = Real(nlayers)  [ BOA output only ]
!  These checks have been implemented inside VLIDORT.

!  1/31/21. Version 2.8.3, Set the NSTOKES2 FOURIER0 flag by hand

      VLIDORT_FixIn%Bool%TS_DO_FOURIER0_NSTOKES2 = .true.

!  Redefine some VLIDORT input variables (override VLIDORT config file)

      VLIDORT_FixIn%Optical%TS_ATMOS_WAVELENGTH = lambdas(w)/1.0D3 !must be in [um]

!  Initialize some input variables not handled by VLIDORT_INPUT_MASTER

      !CALL VLIDORT_BRDF_Sup_Init ( VLIDORT_Sup )
      !CALL VLIDORT_SS_Sup_Init ( VLIDORT_Sup )

!   IMPORTANT - Check compatibility of VSLEAVE and Main VLIDORT inputs
!   ------------------------------------------------------------------

      CALL VLIDORT_VSLEAVE_INPUT_CHECK ( &
        VSLEAVE_Sup_In,             & ! Inputs
        VLIDORT_FixIn,              & ! Inputs
        VLIDORT_ModIn,              & ! Inputs
        VLIDORT_VSLEAVECheck_Status ) ! Outputs

      IF ( VLIDORT_VSLEAVECheck_Status%TS_STATUS_INPUTCHECK .ne. VLIDORT_SUCCESS ) &
        CALL VLIDORT_VSLEAVE_INPUT_CHECK_ERROR ( &
           'V2p8p3_VLIDORT_VSLEAVEcheck_fd_call_Test2.log', VLIDORT_VSLEAVECheck_Status )

!  Finite diff jacobian #1: wrt Pigment Concentration
!  --------------------------------------------------

!  Perturb Pigment Concentration, leave windspeed

      VSLEAVE_Sup_In%SL_CHLORCONC = par(1)*epsfac1
      VSLEAVE_Sup_In%SL_WINDSPEED = par(2)

!  VSLEAVE call - the output will be used for the VLIDORT calculation

      CALL VSLEAVE_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_MAINMASTER failed, perturbation # 1' ; go to 678
      endif

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std only)

      CALL SET_VLIDORT_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out, VLIDORT_FixIn, VLIDORT_ModIn, & !Inputs
        VLIDORT_Sup )                                    !Outputs

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo                        = lambertian_albedo

!  VLIDORT call
!  1/31/21. Version 2.8.3. Need to set the debug_input flag

      do_debug_input = .false.
      CALL VLIDORT_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup, &
        VLIDORT_Out )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save perturbed results

      q = 1
      do v = 1, ngeometries
        I_pert(q,w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
      enddo

!  Finite diff jacobian #2: wrt Windspeed
!  ----------------------------------------

!  Perturb Windspeed, save Pigment Concentration

      VSLEAVE_Sup_In%SL_CHLORCONC = par(1)
      VSLEAVE_Sup_In%SL_WINDSPEED = par(2)*epsfac2

!  VSLEAVE call - the output will be used for the VLIDORT calculation

      CALL VSLEAVE_MAINMASTER ( &
        VSLEAVE_Sup_In,         & ! Inputs
        VSLEAVE_Sup_Out,        & ! Outputs
        VSLEAVE_Sup_OutputStatus) ! Outputs

      if ( VSLEAVE_Sup_OutputStatus%SL_STATUS_OUTPUT .ne. vlidort_success ) THEN
         TRACE = 'VSLEAVE_MAINMASTER failed, perturbation # 2' ; go to 678
      endif

!  Copy VSLEAVE Sup outputs to VLIDORT's VSLEAVE Sup inputs (std only)

      CALL SET_VLIDORT_VSLEAVE_INPUTS ( &
        VSLEAVE_Sup_Out, VLIDORT_FixIn, VLIDORT_ModIn, & !Inputs
        VLIDORT_Sup )                                    !Outputs

!  Set VLIDORT Type structures

      VLIDORT_FixIn%Cont%TS_nlayers                   = nlayers
      VLIDORT_ModIn%MCont%TS_ngreek_moments_input     = ngreek_moments_input
      VLIDORT_FixIn%Chapman%TS_height_grid(0:nlayers) = height_grid(0:nlayers)

      VLIDORT_FixIn%Optical%TS_deltau_vert_input(1:nlayers)             = deltau_vert_input(1:nlayers)
      VLIDORT_ModIn%MOptical%TS_omega_total_input(1:nlayers)            = omega_total_input(1:nlayers)
      VLIDORT_FixIn%Optical%TS_greekmat_total_input(0:nmi,1:nlayers,1)  = greekmat_total_input(0:nmi,1:nlayers,1)
      VLIDORT_FixIn%Optical%TS_lambertian_albedo                        = lambertian_albedo

!  VLIDORT call
!  1/31/21. Version 2.8.3. Need to set the debug_input flag

      do_debug_input = .false.
      CALL VLIDORT_MASTER ( do_debug_input, &
        VLIDORT_FixIn, &
        VLIDORT_ModIn, &
        VLIDORT_Sup, &
        VLIDORT_Out )

!  Exception handling, write-up (optional)

      CALL VLIDORT_WRITE_STATUS ( &
        'V2p8p3_VLIDORT_Execution.log', VLIDORT_ERRUNIT, OPENFILEFLAG, VLIDORT_Out%Status )

!  Save perturbed results

      q = 2
      do v = 1, ngeometries
        I_pert(q,w,v) = VLIDORT_Out%Main%TS_stokes(uta,v,o1,upidx)
      enddo

!  End wavelength loop

   enddo

!  Write results, Test 2
!  ---------------------

   if (do_waterleaving) then
     c2 = 'water'
   else
     c2 = 'clean'
   end if

!  write results to screen
!  -- 5/24/21. revised output.

!   write(*,'(/A/)')'  SCREEN OUTPUT '
!   do w = 1, nlambdas
!      do v = 1, ngeometries
!      write(*,54) lambdas(w),v,&
!         I_base(w,v), ( I_pert(q,w,v), par(q), Jac(q,w,v), par(q)*Jac(q,w,v), &
!                        (I_pert(q,w,v)-I_base(w,v))/eps, q=1,n_totalsurface_wfs )
!     enddo
!   enddo
!54 format(f10.4,i4, es16.8,5es16.8,/,30x,5es16.8)

!  Write Results to file

   open(unit=777,file='vlidort_v_test/results_vsleaveplus_tester_' &
        // trim(c1) // '_' // c2 // '.all',status='replace',action='write')
   write(777,'(/A/A/)')' Wavelength Geom#     Radiance                 Jacobian 1                          Jacobian 2',&
                       '                      Baseline            Analytic  Finite-Diff               Analytic  Finite-Diff'
   do w = 1, nlambdas
      do v = 1, ngeometries
      write(777,55) lambdas(w),v,&
         I_base(w,v), Jac(1,w,v), (I_pert(1,w,v)-I_base(w,v))/eps/par(1), &
                      Jac(2,w,v), (I_pert(2,w,v)-I_base(w,v))/eps/par(2)
     enddo
   enddo
55 format(f10.4,3x,i2,2x,es16.8,2(3x,2es16.8))
   close(777)

!  Finish Test 2

   write(*,*)
   write(*,*) ' TEST 2 finished successfully'

!  continuation point for End of program

89 continue

!  Finish Main

   write(*,*)
   write(*,*) 'Main program finished successfully'
   stop

!mick fix 9/19/2017 - added error finish section
!  Error finish

678   continue

   write(*,*)
   write(*,'(1x,a)') trim(TRACE)
   write(*,*)'Number of error messages from VSLEAVE calculation = ',VSLEAVE_Sup_OutputStatus%SL_NOUTPUTMESSAGES
   do i = 1, VSLEAVE_Sup_OutputStatus%SL_NOUTPUTMESSAGES
      write(*,'(a,i2,a,a)')' Message # ', i, ': ',adjustl(trim(VSLEAVE_Sup_OutputStatus%SL_OUTPUTMESSAGES(i)))
   enddo

   stop
end program VSLEAVEplus_Tester

!

SUBROUTINE RAYLEIGH_FUNCTION &
          ( FORWARD_MAXLAMBDAS, CO2_PPMV_MIXRATIO, FORWARD_NLAMBDAS, FORWARD_LAMBDAS, &
            RAYLEIGH_XSEC, RAYLEIGH_DEPOL )

!  D. Haffner (SSAI) and V. Garcia (DLR). Independently spotted slight bug in coding.
!    -- variable MCO2 is defined in parts per volume (e.g. 0.00036 for 360 ppm).
!    -- OK for Eq.(19) in Bodhaine, but for Eq.(23) the CO2 is in ppv by percent (0.036 for 360 ppm).
!    -- Multiply MCO2 by a factor of 100 in Eq.(23) corrects the error

!  Stand-alone utility routine.
!    -- Rayleigh function (cross-sections and depolarization ratios)
!    -- Bodhaine et. al. (1999) formulae
!    -- CO2 Mixing ratio in PPMV ( e.g. 400).
!    -- Wavelengths in nm (formerly Angstroms)

      IMPLICIT NONE

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

!  wavelength

      INTEGER          FORWARD_MAXLAMBDAS, FORWARD_NLAMBDAS
      DOUBLE PRECISION FORWARD_LAMBDAS ( FORWARD_MAXLAMBDAS )

!  CO2 mixing ratio

      DOUBLE PRECISION CO2_PPMV_MIXRATIO

!  Output arguments
!  ----------------

!  cross-sections and depolarization output

      DOUBLE PRECISION RAYLEIGH_XSEC  ( FORWARD_MAXLAMBDAS )
      DOUBLE PRECISION RAYLEIGH_DEPOL ( FORWARD_MAXLAMBDAS )

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

      INTEGER          W
      DOUBLE PRECISION MASS_DRYAIR
      DOUBLE PRECISION NMOL, PI, CONS
      DOUBLE PRECISION MO2,MN2,MARG,MCO2,MAIR
      DOUBLE PRECISION FO2,FN2,FARG,FCO2,FAIR
      DOUBLE PRECISION LAMBDA_C,LAMBDA_M,LPM2,LP2
      DOUBLE PRECISION N300M1,NCO2M1,NCO2,XCO2
      DOUBLE PRECISION NCO2SQ, NSQM1,NSQP2,TERM

!  data statements and parameters
!  ------------------------------

      DATA MO2  / 20.946D0 /
      DATA MN2  / 78.084D0 /
      DATA MARG / 0.934D0 /

      DOUBLE PRECISION, PARAMETER :: S0_A = 15.0556D0
      DOUBLE PRECISION, PARAMETER :: S0_B = 28.9595D0

      DOUBLE PRECISION, PARAMETER :: S1_A = 8060.51D0
      DOUBLE PRECISION, PARAMETER :: S1_B = 2.48099D+06
      DOUBLE PRECISION, PARAMETER :: S1_C = 132.274D0
      DOUBLE PRECISION, PARAMETER :: S1_D = 1.74557D+04
      DOUBLE PRECISION, PARAMETER :: S1_E = 39.32957D0

      DOUBLE PRECISION, PARAMETER :: S2_A = 0.54D0

      DOUBLE PRECISION, PARAMETER :: S3_A = 1.034D0
      DOUBLE PRECISION, PARAMETER :: S3_B = 3.17D-04
      DOUBLE PRECISION, PARAMETER :: S3_C = 1.096D0
      DOUBLE PRECISION, PARAMETER :: S3_D = 1.385D-03
      DOUBLE PRECISION, PARAMETER :: S3_E = 1.448D-04

!  Start of code
!  -------------

!  constants

      NMOL = 2.546899D19
      PI   = DATAN(1.0D0)*4.0D0
      CONS = 24.0D0 * PI * PI * PI

!  convert co2. MCO2 is in Parts per Volume

      MCO2 = 1.0D-06 * CO2_PPMV_MIXRATIO

!  mass of dry air: Eq.(17) of BWDS

      MASS_DRYAIR = S0_A * MCO2 + S0_B

!  start loop

      DO W = 1, FORWARD_NLAMBDAS

!  wavelength in micrometers

      LAMBDA_M = 1.0D-03 * FORWARD_LAMBDAS(W)
      LAMBDA_C = 1.0D-07 * FORWARD_LAMBDAS(W)
      LPM2     = 1.0D0 / LAMBDA_M / LAMBDA_M

!  step 1: Eq.(18) of BWDS

      N300M1 = S1_A + ( S1_B / ( S1_C - LPM2 ) ) + &
                      ( S1_D / ( S1_E - LPM2 ) )
      N300M1 = N300M1 * 1.0D-08

!  step 2: Eq.(19) of BWDS. Uses MCO2, Parts per Volume.

      NCO2M1 = N300M1 * ( 1.0D0 + S2_A * ( MCO2  - 0.0003D0 ) )
      NCO2   = NCO2M1 + 1
      NCO2SQ = NCO2 * NCO2

!  step 3: Eqs. (5&6) of BWDS (Bates' results)

      FN2  = S3_A + S3_B * LPM2
      FO2  = S3_C + S3_D * LPM2 + S3_E * LPM2 * LPM2

!  step 4: Eq.(23) of BWDS. Uses Parts per Volume in percent, i.e. XCO2 = 100 * MCO2
!     ---> King factor and depolarization ratio

      FARG = 1.0D0
      FCO2 = 1.15D0

!mick note 4/20/2021 - old lines turned back on to compare with VLIDORT

!  Old lines
      MAIR = MN2 + MO2 + MARG + MCO2
      FAIR = MN2*FN2 + MO2*FO2 + MARG*FARG + MCO2*FCO2

!  New lines. Use XCO2 instead of MCO2
!      XCO2 = 100.0_fpk * MCO2
!      MAIR = MN2 + MO2 + MARG + XCO2
!      FAIR = MN2*FN2 + MO2*FO2 + MARG*FARG + XCO2 *FCO2

      FAIR = FAIR / MAIR
      RAYLEIGH_DEPOL(W) = 6.0D0*(FAIR-1.0D0)/(3.0D0+7.0D0*FAIR)

!  step 5: Eq.(22) of BWDS
!     ---> Cross section

      LP2  = LAMBDA_C * LAMBDA_C
      NSQM1 = NCO2SQ - 1.0D0
      NSQP2 = NCO2SQ + 2.0D0
      TERM = NSQM1 / LP2 / NMOL / NSQP2
      RAYLEIGH_XSEC(W) =  CONS * TERM * TERM * FAIR

!  end loop

      ENDDO

!  finish
!  ------

      RETURN
END SUBROUTINE RAYLEIGH_FUNCTION

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps