text:
*DECK J4SAVE
      FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
C***BEGIN PROLOGUE  J4SAVE
C***SUBSIDIARY
C***PURPOSE  Save or recall global variables needed by error
C            handling routines.
C***LIBRARY   SLATEC (XERROR)
C***TYPE      INTEGER (J4SAVE-I)
C***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        J4SAVE saves and recalls several global variables needed
C        by the library error handling routines.
C
C     Description of Parameters
C      --Input--
C        IWHICH - Index of item desired.
C                = 1 Refers to current error number.
C                = 2 Refers to current error control flag.
C                = 3 Refers to current unit number to which error
C                    messages are to be sent.  (0 means use standard.)
C                = 4 Refers to the maximum number of times any
C                     message is to be printed (as set by XERMAX).
C                = 5 Refers to the total number of units to which
C                     each error message is to be written.
C                = 6 Refers to the 2nd unit for error messages
C                = 7 Refers to the 3rd unit for error messages
C                = 8 Refers to the 4th unit for error messages
C                = 9 Refers to the 5th unit for error messages
C        IVALUE - The value to be set for the IWHICH-th parameter,
C                 if ISET is .TRUE. .
C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
C                 given the value, IVALUE.  If ISET=.FALSE., the
C                 IWHICH-th parameter will be unchanged, and IVALUE
C                 is a dummy parameter.
C      --Output--
C        The (old) value of the IWHICH-th parameter will be returned
C        in the function value, J4SAVE.
C
C***SEE ALSO  XERMSG
C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
C                 Error-handling Package, SAND82-0800, Sandia
C                 Laboratories, 1982.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900205  Minor modifications to prologue.  (WRB)
C   900402  Added TYPE section.  (WRB)
C   910411  Added KEYWORDS section.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  J4SAVE
      LOGICAL ISET
      INTEGER IPARAM(9)
      SAVE IPARAM
      DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,6,10/
      DATA IPARAM(5)/1/
      DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
C***FIRST EXECUTABLE STATEMENT  J4SAVE
      J4SAVE = IPARAM(IWHICH)
      IF (ISET) IPARAM(IWHICH) = IVALUE
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


use,intrinsic:: iso_fortran_env, only: sp=>real32, dp=>real64, qp=>real128
implicit none

integer :: ps,pd,pq
real(sp) :: eps32,  s
real(dp) :: eps64,  d
real(qp) :: eps128, q

! ----- single prec -------
ps=0
s=1.

do while (s /= s + 1.)
  s = s * 2.
  ps = ps + 1
end do

eps32 = 2.**(-(ps-1.))
print *, 'bits of precision for 32-bit float', ps
print *, 'machine epsilon: 32-bit float', eps32
! ----- double prec -------
pd=0
d=1._dp

do while (d /= d + 1._dp)
  d = d * 2._dp
  pd = pd + 1
end do

eps64 = 2._dp**(-(pd-1._dp))
print *, 'bits of precision for 64-bit float', pd
print *, 'machine epsilon: 64-bit float', eps64
! ----- quad prec -------
pq=0
q=1._qp

do while (q /= q + 1._qp)
  q = q * 2._qp
  pq = pq + 1
end do

eps128 = 2._qp**(-(pq-1._qp))
print *, 'bits of precision for 128-bit float', pq
print *, 'machine epsilon: 128-bit float', eps128

end program


text:
C  =====================================================================
C  pgm: SHPCOD .. Get parameter code from message
C
C  use:     CALL SHPCOD(KHAR,KHPOS)
C
C  i/o: KHAR ...... last buffer char obtained - CHAR*1
C  i/o: KHPOS ..... last char loc: 2=eol,1=err-eol,0=eof,neg=err - INT
C
C  rqd: SHGETK,SH2BLA,SHERR,SHSAVP
C
C  cmt: Routine will NOT process if current char is "end-of-line".
C  cmt: "KHAR" must be guaranteed not blank if "KHPOS" is less than 2.
C  =====================================================================
      SUBROUTINE SHPCOD(KHAR,KHPOS)

      EXTERNAL       SHGETK,SH2BLA,SHERR,SHSAVP

      CHARACTER*8    PARCOD
      CHARACTER*1    KHAR
      INTEGER        KHPOS,NOC,NN,III
C
C    ================================= RCS keyword statements ==========
      CHARACTER*68     RCSKW1,RCSKW2
      DATA             RCSKW1,RCSKW2 /                                 '
     .$Source: /fs/hseb/ob5/rfc/ofs/src/shefpars_driv/RCS/shpcod.f,v $
     . $',                                                             '
     .$Id: shpcod.f,v 1.3 1998/04/07 19:14:39 page Exp $
     . $' /
C    ===================================================================
C

        CALL SH2BLA(KHAR,KHPOS,NOC)

C                   If current char position is in a message line

        IF (KHPOS .GT. 2) THEN

          CALL SHSAVP('G',III,PARCOD)

C                   Loop while char is not blank nor is output str full

          NN = 0
   30     IF(KHPOS.LE.2.OR.NN.GE.7.OR.KHAR.EQ.' '.OR.KHAR.EQ.'/')GOTO 40

            IF (KHAR.GE.'A' .AND. KHAR.LE.'Z') THEN
              NN = NN+1
              PARCOD(NN:NN) = KHAR
              CALL SHGETK(KHAR,KHPOS)
            ELSEIF (NN.EQ.3 .OR. NN.EQ.4 .OR. NN.EQ.6) THEN
              IF (KHAR.GE.'1' .AND. KHAR.LE.'9') THEN
                NN = NN+1
                PARCOD(NN:NN) = KHAR
                CALL SHGETK(KHAR,KHPOS)
              ELSE
                CALL SHERR('E',29,KHPOS,KHAR)
              ENDIF
            ELSE
              CALL SHERR('E',29,KHPOS,KHAR)
            ENDIF

            GOTO 30
   40     CONTINUE

          CALL SHSAVP('P',III,PARCOD)

          CALL SH2BLA(KHAR,KHPOS,NOC)

          IF (KHPOS.GT.2 .AND. KHAR.NE.'/') THEN
            IF (NOC .EQ. 0) CALL SHERR('E',54,KHPOS,KHAR)
          ENDIF

          IF (KHPOS.NE.1 .AND. PARCOD(2:2).EQ.'        ') THEN
            CALL SHERR('E',64,KHPOS,KHAR)
          ENDIF

        ENDIF

      RETURN
      END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


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

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

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

module lidort_geometry

!  Parameter types

   USE LIDORT_PARS, only : fpk

!private
public

! ###############################################################
! #                                                             #
! # Subroutines in this Module                                  #
! #                                                             #
! #                                                             #
! #            LIDORT_CHAPMAN   (called by LIDORT_MASTER)       #
! #              BEAM_GEOMETRY_PREPARE                          #
! #                                                             #
! #            OUTGOING_SPHERGEOM_FINE_UP                       #
! #            OUTGOING_SPHERGEOM_FINE_DN                       #
! #            MULTI_OUTGOING_ADJUSTGEOM                        #
! #            OBSGEOM_OUTGOING_ADJUSTGEOM                      #
! #            LOSONLY_OUTGOING_ADJUSTGEOM                      #
! #                                                             #
! ###############################################################

! ###############################################################
! #                                                             #
! #  Revisions in Version 3.5                                   #
! #                                                             #
! #     07 October 2010.                                        #
! #       Original code had only one geometry routine for the   #
! #       Outgoing spherical geometry computation. This was     #
! #       only correct for the upwelling case, and actually     #
! #       incorrect (in a small way) for the (rarely-used)      #
! #       downwelling case. Now, there are separate outgoing    #
! #       sphergeom routine for "_up" and "dn" and these are    #
! #       called separately in the various SSCORR_OUTGOING      #
! #       single scatter routines.                              #
! #                                                             #
! ###############################################################

contains

SUBROUTINE LIDORT_CHAPMAN                                 &
          ( DO_PLANE_PARALLEL, DO_REFRACTIVE_GEOMETRY,    & !  Input
            NLAYERS, NBEAMS, FINEGRID, BEAM_SZAS,         & !  Input
            EARTH_RADIUS, RFINDEX_PARAMETER,              & !  Input
            HEIGHT_GRID, PRESSURE_GRID, TEMPERATURE_GRID, & !  Input
            CHAPMAN_FACTORS, SZA_LOCAL_INPUT,             & !  Output
            FAIL, MESSAGE, TRACE )                          !  Output

!  The following options apply:

!   1. If the plane-parallel flag is on, no further inputs are required

!   2. If the plane-parallel flag is off, then Pseudo-spherical:
!       (a) Straight line geometry, must specify
!               Earth_radius, height grid
!       (b) Refractive geometry, must specify
!               Earth_radius, height grid
!               pressure grid, temperature grid

!  The logic will be checked before the module is called.

!  Newly programmed by R. Spurr, RT SOLUTIONS Inc. 5/5/05.

!    Based round a call to a pure geometry module which returns slant
!    path distances which was adapted for use in the Radiant model by
!    R. Spurr during an OCO L2 intensive April 24-29, 2005.

!  module, dimensions and numbers

      USE LIDORT_pars, only : MAXBEAMS, MAXLAYERS

!  implicit none

      IMPLICIT NONE

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

!  flag for plane parallel case

      LOGICAL         :: DO_PLANE_PARALLEL

!  flag for refractive geometry

      LOGICAL         :: DO_REFRACTIVE_GEOMETRY

!  Number of layers

      INTEGER         :: NLAYERS

!  number of solar beams to be processed

      INTEGER         :: NBEAMS

!  number of fine layers within coarse layers

      INTEGER         :: FINEGRID(MAXLAYERS)

!  TOA solar zenith angles

      REAL(fpk)       :: BEAM_SZAS ( MAXBEAMS )

!  Earth radius (km)

      REAL(fpk)       :: EARTH_RADIUS

!  Refractive index parametaer (Born-Wolf approximation)

      REAL(fpk)       :: RFINDEX_PARAMETER

!  Coarse grids of heights, pressures and temperatures

      REAL(fpk)       :: HEIGHT_GRID     (0:MAXLAYERS)
      REAL(fpk)       :: PRESSURE_GRID   (0:MAXLAYERS)
      REAL(fpk)       :: TEMPERATURE_GRID(0:MAXLAYERS)

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

!  Chapman factors

      REAL(fpk)       :: CHAPMAN_FACTORS ( MAXLAYERS, MAXLAYERS, MAXBEAMS )

!  solar zenith angles at nadir

      REAL(fpk)       :: SZA_LOCAL_INPUT(0:MAXLAYERS,MAXBEAMS)

!  output status

      LOGICAL         :: FAIL
      CHARACTER*(*)   :: MESSAGE, TRACE

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

!  number of iterations (refractive case only)
!      This is debug output

      INTEGER         :: ITERSAVE(MAXLAYERS)

!  other local variables

      INTEGER      :: IB
      REAL(fpk)    :: SUN0

!  get spherical optical depths
!  ----------------------------

!  start beam loop

      DO IB = 1, NBEAMS

        SUN0 = BEAM_SZAS(IB)

        CALL BEAM_GEOMETRY_PREPARE                         &
          ( MAXBEAMS, MAXLAYERS, IB, NLAYERS, FINEGRID,    & ! Input
            SUN0, EARTH_RADIUS, RFINDEX_PARAMETER,         & ! Input
            DO_PLANE_PARALLEL, DO_REFRACTIVE_GEOMETRY,     & ! Input
            HEIGHT_GRID, PRESSURE_GRID, TEMPERATURE_GRID,  & ! Input
            CHAPMAN_FACTORS, SZA_LOCAL_INPUT,             & ! Input
            ITERSAVE, FAIL, MESSAGE )                        ! Output

!  return if failed

        IF ( FAIL ) THEN
          TRACE = ' BEAM_GEOMETRY_PREPARE Geometry failure'
          RETURN
        ENDIF

!  end beam loop

      ENDDO

!  Finish

      RETURN
END SUBROUTINE LIDORT_CHAPMAN

!

SUBROUTINE BEAM_GEOMETRY_PREPARE                           &
          ( MAXBEAMS, MAXLAYERS, IBEAM, NLAYERS, FINEGRID, & ! Input
            SZA_GEOM_TRUE, REARTH, RFINDEX_PARAMETER,      & ! Input
            DO_PLANE_PARALLEL, DO_REFRACTIVE_GEOMETRY,     & ! Input
            HEIGHTS, PRESSURES, TEMPERATURES,              & ! Input
            CHAPMAN_FACTORS, SZA_LEVEL_OUTPUT,             & ! Input
            ITERSAVE, FAIL, MESSAGE )                        ! Output

!  Module, dimensions and numbers

      USE LIDORT_pars, only : ZERO, ONE

!  Implicit none

      IMPLICIT NONE

!  Generate path CHAPMAN_FACTORS and SZA angles SZA_LEVEL_OUTPUT
!  for a curved ray-traced beam through a multilayer atmosphere.

!  Coarse layering is input to the module. Values of Z, P, T  are
!  given at the layer boundaries, with the first value (index 0) at TOA.

!  The refractive geometry is assumed to start at the TOA level.

!  We also require the earth radius and the refractive index parameter
!   (For the Born-Wolf approximation)

!  There is no refraction if the flag DO_REFRACTIVE_GEOMETRY is not set.
!  In this case we do not require pressure and temperature information.
!  The calculation will then be for geometric rays.

!  The plane parallel Flag and the refractive geometry flag should not
!  both be true - this should be checked outside.

!  In the refracting case, fine-gridding of pressure and temperature is
!  done internally, temperature is interpolated linearly with height,
!  and pressure log-linearly. The refraction uses Snell's law rule.
!  Finelayer gridding assumes equidistant heights within coarse layers
!  but the number of fine layers can be varied

!  Output is specified at coarse layer boundaries

!  Module is stand-alone.

!  Reprogrammed for the OCO L2 algorithm
!   R. Spurr, RT Solutions, Inc.   April 27, 2005

!  Intended use in LIDORT and Radiant RT models.

!  Input arguments
!  ===============

!  input dimensioning

      INTEGER      :: MAXLAYERS, MAXBEAMS

!  Beam index

      INTEGER      :: IBEAM

!  number of coarse layers

      INTEGER      :: NLAYERS

!  number of fine layers within coarse layers

      INTEGER      :: FINEGRID(MAXLAYERS)

!  True solar zenith angle (degrees)

      REAL(fpk)    :: SZA_GEOM_TRUE

!  Earth radius (km)

      REAL(fpk)    :: REARTH

!  Refractive index parametaer (Born-Wolf approximation)

      REAL(fpk)    :: RFINDEX_PARAMETER

!  flag for plane parallel case

      LOGICAL      :: DO_PLANE_PARALLEL

!  flag for refractive geometry

      LOGICAL      :: DO_REFRACTIVE_GEOMETRY

!  Coarse grids of heights, pressures and temperatures

      REAL(fpk)    :: HEIGHTS     (0:MAXLAYERS)
      REAL(fpk)    :: PRESSURES   (0:MAXLAYERS)
      REAL(fpk)    :: TEMPERATURES(0:MAXLAYERS)

!  Output arguments
!  ================

!  Path segments distances (km)

      REAL(fpk)      :: CHAPMAN_FACTORS(MAXLAYERS,MAXLAYERS,MAXBEAMS)

!  solar zenith angles at nadir

      REAL(fpk)      :: SZA_LEVEL_OUTPUT(0:MAXLAYERS,MAXBEAMS)

!  number of iterations (refractive case only)
!   This is debug output

      INTEGER         :: ITERSAVE(MAXLAYERS)

!  output status

      LOGICAL         :: FAIL
      CHARACTER*(*)   :: MESSAGE

!  Local variables
!  ===============

!  local dimensioning

      INTEGER, PARAMETER ::  LOCAL_MAXLAYERS     = 200
      INTEGER, PARAMETER ::  LOCAL_MAXFINELAYERS = 10

!  fine layer gridding for refraction

      REAL(fpk)    :: ZRFINE(LOCAL_MAXLAYERS,0:LOCAL_MAXFINELAYERS)
      REAL(fpk)    :: PRFINE(LOCAL_MAXLAYERS,0:LOCAL_MAXFINELAYERS)
      REAL(fpk)    :: TRFINE(LOCAL_MAXLAYERS,0:LOCAL_MAXFINELAYERS)

!  local height arrays

      REAL(fpk)    :: H(0:LOCAL_MAXLAYERS)
      REAL(fpk)    :: DELZ(LOCAL_MAXLAYERS)

!  help variables

      INTEGER      :: N, J, NRFINE, K, ITER, MAXF,IB
      LOGICAL      :: LOOP
      REAL(fpk)    :: GM_TOA, TH_TOA, MU_TOA, MU_NEXT
      REAL(fpk)    :: Z1, Z0, Z, T1, T0, T, P1, P0, Q1, Q0, Q
      REAL(fpk)    :: FU, FL, DEG_TO_RAD

      REAL(fpk)    :: LAYER_DIST, MU_PREV, STH1, SINTH1, STH2, SINTH2, LOCAL_SUBTHICK, &
                      PHI, PHI_0, PHI_CUM, SINPHI, DELPHI, REFRAC, RATIO,              &
                      RE_LOWER, RE_UPPER, DIST, STH2D, SINTH2D, SNELL

!  Standard temperature (K) and pressure (mbar).

      REAL(fpk), PARAMETER  :: T_STANDARD = 273.16D0
      REAL(fpk), PARAMETER  :: P_STANDARD = 1013.25D0
      REAL(fpk), PARAMETER  :: STP_RATIO = T_STANDARD / P_STANDARD

!  Loschmidt's number (particles/cm2/km).

      REAL(fpk), PARAMETER  :: RHO_STANDARD = 2.68675D+24

!  Some setup operations
!  =====================

!  initialise output

      IB = IBEAM
      SZA_LEVEL_OUTPUT(0,IB) = 0.0D0
      DO N = 1, NLAYERS
        SZA_LEVEL_OUTPUT(N,IB) = 0.0D0
        DO K = 1, NLAYERS
          CHAPMAN_FACTORS(N,K,IB) = 0.0D0
        ENDDO
      ENDDO

      FAIL    = .FALSE.
      MESSAGE = ' '

!  check local dimensioning

      IF ( LOCAL_MAXLAYERS .LT. NLAYERS ) THEN
        MESSAGE = 'local coarse layer dimensioning insufficient'
        FAIL = .TRUE.
        RETURN
      ENDIF

!  Check fine layers do not exceed local dimensions assigned

      IF ( DO_REFRACTIVE_GEOMETRY ) THEN
       MAXF = 0
       DO N = 1, NLAYERS
         MAXF = MAX(MAXF,FINEGRID(N))
       ENDDO
       IF ( LOCAL_MAXFINELAYERS .LT. MAXF ) THEN
         MESSAGE = 'local fine layer dimensioning insufficient'
         FAIL = .TRUE.
         RETURN
       ENDIF
      ENDIF

!  earth radii and heights differences

      DO N = 0, NLAYERS
        H(N) = HEIGHTS(N) + REARTH
      ENDDO

      DO N = 1, NLAYERS
        DELZ(N) = HEIGHTS(N-1)-HEIGHTS(N)
      ENDDO

!  TOA values

      SZA_LEVEL_OUTPUT(0,IB) = SZA_GEOM_TRUE
      DEG_TO_RAD = DATAN(1.0D0) / 45.0D0
      TH_TOA = SZA_GEOM_TRUE * DEG_TO_RAD
      MU_TOA = DCOS(TH_TOA)
      GM_TOA = DSQRT ( 1.0D0 - MU_TOA * MU_TOA )

!  initialize

      STH2D  = 0.0D0

!  derive the fine values

      IF ( DO_REFRACTIVE_GEOMETRY ) THEN

        Z0 = HEIGHTS(0)
        P0 = PRESSURES(0)
        T0 = TEMPERATURES(0)
        Q0 = DLOG(P0)
        DO N = 1, NLAYERS
          NRFINE = FINEGRID(N)
          LOCAL_SUBTHICK = DELZ(N) / DBLE(NRFINE)
          P1 = PRESSURES(N)
          Z1 = HEIGHTS(N)
          T1 = TEMPERATURES(N)
          Q1 = DLOG(P1)
          ZRFINE(N,0) = Z0
          PRFINE(N,0) = P0
          TRFINE(N,0) = T0
          DO J = 1, NRFINE - 1
            Z  = Z0 - DBLE(J)*LOCAL_SUBTHICK
            FL = ( Z0 - Z ) / DELZ(N)
            FU = 1.0d0 - FL
            Q  = FL * Q1 + FU * Q0
            T  = FL * T0 + FU * T1
            PRFINE(N,J) = DEXP (  Q )
            TRFINE(N,J) = T
            ZRFINE(N,J) = Z
          ENDDO
          PRFINE(N,NRFINE) = P1
          TRFINE(N,NRFINE) = T1
          ZRFINE(N,NRFINE) = Z1
!              write(*,'(i3,11F10.4)')N,(PRFINE(N,J),J=0,NRFINE)
          Z0 = Z1
          P0 = P1
          T0 = T1
          Q0 = Q1
        ENDDO
      ENDIF

!  plane-parallel case
!  ===================

      IF ( DO_PLANE_PARALLEL ) THEN
        DO N = 1, NLAYERS
          SZA_LEVEL_OUTPUT(N,IB) = SZA_GEOM_TRUE
          DO K = 1, N
            CHAPMAN_FACTORS(N,K,IB) = 1.0D0 / MU_TOA
          ENDDO
        ENDDO
        RETURN
      ENDIF

!  Refractive Geometry case
!  ========================

      IF ( DO_REFRACTIVE_GEOMETRY ) THEN

!  Rob fix 9/9/14. Sun zero case

        IF ( GM_TOA .eq. zero ) then
          DO N = 1, NLAYERS
            CHAPMAN_FACTORS(N,1:N,IB) = one
          ENDDO
          RETURN
        ENDIF

!  Non-zero case....(resume code)

!  Start value of SZA cosine

        MU_PREV = MU_TOA

!  start layer loop

        DO N = 1, NLAYERS

!  start values

          SINTH1 = GM_TOA * H(N) / H(0)
          STH1 = DASIN(SINTH1)
          PHI_0 = TH_TOA - STH1
          NRFINE = FINEGRID(N)

!  iteration loop

          ITER = 0
          LOOP = .TRUE.
          DO WHILE (LOOP.AND.ITER.LT.100)
            ITER = ITER + 1
            PHI_CUM = 0.0D0
            RE_UPPER = ZRFINE(1,0) + REARTH
            RATIO  = PRFINE(1,0) * STP_RATIO / TRFINE(1,0)
            REFRAC = 1.0D0 + RFINDEX_PARAMETER * RATIO
            SNELL = REFRAC * RE_UPPER * SINTH1
            DO K = 1, N
              LAYER_DIST = 0.0D0
              LOCAL_SUBTHICK = DELZ(K) / DBLE(NRFINE)
              DO J = 0, NRFINE - 1
                RATIO  = PRFINE(K,J) * STP_RATIO / TRFINE(K,J)
                REFRAC = 1.0D0 + RFINDEX_PARAMETER * RATIO
                RE_LOWER = RE_UPPER - LOCAL_SUBTHICK
                SINTH2 = SNELL/ (REFRAC * RE_UPPER )
                IF ( SINTH2.GT.1.0D0 ) SINTH2 = 1.0D0
                STH2 = DASIN(SINTH2)
                SINTH2D = RE_UPPER * SINTH2 / RE_LOWER
                IF ( SINTH2D .GT. 1.0D0 ) THEN
                  MESSAGE = 'refraction yields angles > 90 some levels'
                  FAIL = .TRUE.
                  RETURN
                ENDIF
                STH2D = DASIN(SINTH2D)
                PHI = STH2D - STH2
                SINPHI = DSIN(PHI)
                PHI_CUM = PHI_CUM + PHI
                DIST = RE_UPPER * SINPHI / SINTH2D
                LAYER_DIST = LAYER_DIST +  DIST
                RE_UPPER = RE_LOWER
              ENDDO
              CHAPMAN_FACTORS(N,K,IB) = LAYER_DIST / DELZ(K)
            ENDDO

!  examine convergence

            DELPHI = PHI_0 - PHI_CUM
            LOOP = (DABS(DELPHI/PHI_CUM).GT.1.0D-4)

!  Fudge factors to speed up the iteration

            IF ( SZA_GEOM_TRUE .GT. 88.7D0 ) THEN
              STH1 = STH1 + 0.1 * DELPHI
              PHI_0 = TH_TOA - STH1
            ELSE IF ( SZA_GEOM_TRUE .LT. 80.0D0 ) THEN
              PHI_0 = PHI_CUM
              STH1 = TH_TOA - PHI_0
            ELSE
              STH1 = STH1 + 0.3 * DELPHI
              PHI_0 = TH_TOA - STH1
            ENDIF
            SINTH1 = DSIN(STH1)

          ENDDO

!  failure

          IF ( LOOP ) THEN
            MESSAGE = 'refractive iteration not converged'
            FAIL = .TRUE.
            RETURN
          ENDIF

!  Update and save angle output

          MU_NEXT = DCOS(STH2D)
          MU_PREV = MU_NEXT
          SZA_LEVEL_OUTPUT(N,IB) = DACOS(MU_NEXT) / DEG_TO_RAD
          ITERSAVE(N) = ITER

        ENDDO

!  Straight line geometry
!  ======================

      ELSE

!  Rob fix 9/9/14. Sun zero case

        IF ( GM_TOA .eq. zero ) then
          DO N = 1, NLAYERS
            CHAPMAN_FACTORS(N,1:N,IB) = one
          ENDDO
          RETURN
        ENDIF

!  Non-zero case....(resume code)

        DO N = 1, NLAYERS

!  start values

          SINTH1 = GM_TOA * H(N) / H(0)
          STH1   = DASIN(SINTH1)
          RE_UPPER = H(0)

!  solar zenith angles are all the same = input value

          SZA_LEVEL_OUTPUT(N,IB) = SZA_GEOM_TRUE

! loop over layers K from 1 to layer N

          DO K = 1, N

!  sine-rule; PHI = earth-centered angle

            RE_LOWER = RE_UPPER - DELZ(K)
            SINTH2 = RE_UPPER * SINTH1 / RE_LOWER
            STH2   = DASIN(SINTH2)
            PHI    = STH2 - STH1
            SINPHI = DSIN(PHI)
            DIST = RE_UPPER * SINPHI / SINTH2
            CHAPMAN_FACTORS(N,K,IB) = DIST / DELZ(K)

!  re-set

            RE_UPPER = RE_LOWER
            SINTH1 = SINTH2
            STH1   = STH2

          ENDDO

!  finish main layer loop

        ENDDO

!  Finish

      ENDIF

!  end of routine

      RETURN
END SUBROUTINE BEAM_GEOMETRY_PREPARE

!

SUBROUTINE outgoing_sphergeom_fine_up                              &
       ( maxlayers, maxfine, maxpartials,                          & ! Input
         do_fine, do_partials, nlayers, nfine,                     & ! Input
         n_partials, partials_idx,                                 & ! Input
         heights, heights_p, eradius,                              & ! Input
         alpha_boa, theta_boa, phi_boa,                            & ! Input
         sunpaths,      radii,      ntraverse,      alpha_all,     & ! Output
         sunpaths_fine, ntraverse_fine, alpha_fine,                & ! Output
         sunpaths_p,    ntraverse_p,    alpha_p,                   & ! Output
         partials_fineidx, cosscat_up, fail, message )               ! Output

!  Completely stand-alone geometry routine for the outgoing correction
!     This is applicable to the Upwelling path geometry

!    starting inputs are the BOA values of SZA, VZA and PHI
!    need also the height grids, earth radius and control

!  This routine has the fine gridding treatment

!  Version 3.3. September 2007, Partial layer geometries added
!  Version 3.5. Revision to separate upwelling/downwelling geometry routines

      implicit none

!  inputs

      integer       :: maxlayers, maxfine, maxpartials
      integer       :: nlayers, nfine
      logical       :: do_fine, do_partials
      integer       :: n_partials
      integer       :: partials_idx    (maxpartials)
      real(fpk)     :: eradius, heights (0:maxlayers)
      real(fpk)     :: heights_p(maxpartials)
      real(fpk)     :: alpha_boa, theta_boa, phi_boa

!  main outputs (geometry)

      integer       :: ntraverse   (0:maxlayers)
      real(fpk)     :: sunpaths   (0:maxlayers,maxlayers)
      real(fpk)     :: radii      (0:maxlayers)
      real(fpk)     :: alpha_all  (0:maxlayers)

!  Fine level output (geometry)

      integer       :: ntraverse_fine(maxlayers,maxfine)
      real(fpk)     :: sunpaths_fine (maxlayers,maxlayers,maxfine)
      real(fpk)     :: alpha_fine    (maxlayers,maxfine)

!  Partial layer output (geometry)

      integer       :: ntraverse_p(maxpartials)
      real(fpk)     :: sunpaths_p (maxpartials,maxlayers)
      real(fpk)     :: alpha_p    (maxpartials)

!  scattering angles

      real(fpk)     :: cosscat_up (0:maxlayers)

!   Local arrays

      real(fpk)     :: radii_fine    (maxlayers,maxfine)
      real(fpk)     :: radii_p    (maxpartials)
      real(fpk)     :: lospaths(maxlayers)

      real(fpk)     :: lospaths_p_up(maxpartials)
      integer       :: partials_fineidx(maxpartials)

      real(fpk)     :: theta_all  (0:maxlayers)
      real(fpk)     :: phi_all    (0:maxlayers)

!  Status output

      logical       :: fail
      character*(*) :: message

!  Local

      logical       :: direct_sun
      integer       :: n, j, k, krad, n1, ut, np
      real(fpk)     :: deg_to_rad, ex, ey, ez, px, py, pz
      real(fpk)     :: salpha_boa, calpha_boa, sphi_boa
      real(fpk)     :: stheta_boa, ctheta_boa, cphi_boa
      real(fpk)     :: ksi, cksi, sksi, xicum, tangr, fac
      real(fpk)     :: ctheta, stheta, calpha, salpha, cphi
      real(fpk)     :: b, sth0, th0, ks1, sth1, th1, theta_p

!  Local arrays associated with fine grid output

      integer, parameter :: maxlocalfine = 20
      logical            :: direct_sunf(maxlocalfine)
      real(fpk)          :: difz, dfine1, saf, xicum0, path
      real(fpk)          :: thetaf(maxlocalfine), xicumf, difa
      real(fpk)          :: cthetaf(maxlocalfine)
      real(fpk)          :: sthetaf(maxlocalfine)
      real(fpk)          :: ksif(maxlocalfine)

!  Initialise output

      fail = .false.
      message = ' '

!  check range of inputs

      if ( alpha_boa.ge.90.0d0.or.alpha_boa.lt.0.0d0 ) then
        message = 'boa LOS angle outside range [0,90]); Check it!'
        fail    = .true.
        return
      endif
      if ( phi_boa.lt.0.0d0 )   phi_boa = - phi_boa
      if ( phi_boa.gt.180.0d0 ) phi_boa = 360.0d0 - phi_boa
      if ( theta_boa.ge.90.0d0.or.theta_boa.lt.0.0d0 ) then
        message = 'boa SZA angle outside range [0,90]); Check it!'
        fail    = .true.
        return
      endif
      if ( do_fine ) then
       if ( nfine.gt.maxlocalfine ) then
         message = 'local finelayer dimensioning insufficient'
         fail    = .true.
         return
       endif
      endif

!  zero the sun paths
!  Initialize number of layers traversed  (nominal conditions)

      do n = 0, nlayers
        ntraverse(n) = n
        do k = 1, nlayers
         sunpaths(n,k) = 0.0d0
        enddo
      enddo

!  Zero the partial paths, initialize the traverse number (nominal)

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)
          ntraverse_p(ut) = np
          do k = 1, np
            sunpaths_p(ut,k) = 0.0d0
          enddo
        enddo
      endif

!  Zero the fine data paths

      if ( do_fine ) then
       dfine1 = dble(nfine) + 1
       do n = 1, nlayers
        do j = 1, nfine
         ntraverse_fine(n,j) = n
         do k = 1, nlayers
          sunpaths_fine(n,k,j) = 0.0d0
         enddo
        enddo
       enddo
      endif

!  start at BOA

      deg_to_rad = dacos(-1.0d0) / 180.0d0
      alpha_all(nlayers) = alpha_boa * deg_to_rad
      theta_all(nlayers) = theta_boa * deg_to_rad
      phi_all(nlayers)   = phi_boa   * deg_to_rad

!  Cosine of scattering angle at boa

      salpha_boa = dsin(alpha_all(nlayers))
      calpha_boa = dcos(alpha_all(nlayers))
      stheta_boa = dsin(theta_all(nlayers))
      ctheta_boa = dcos(theta_all(nlayers))
      cphi_boa   = dcos(phi_all(nlayers))
      sphi_boa   = dsin(phi_all(nlayers))

      cosscat_up (nlayers) = - calpha_boa * ctheta_boa + &
                               salpha_boa * stheta_boa * cphi_boa

!  Radii
!  -----

!  layer levels

      do n = 0, nlayers
        radii(n) = eradius + heights(n)
      enddo

!  Fine levels

      if ( do_fine ) then
        do n = 1, nlayers
          difz = (radii(n-1)-radii(n))/dfine1
          do j = 1, nfine
            radii_fine(n,j) = radii(n) + difz * dble(j)
          enddo
        enddo
      endif

!  Partial levels
!   Find the fine level just below partial point = (0,1,2,...nfine)

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)
          radii_p(ut) = eradius + heights_p(ut)
          j = 1
          do while (radii_p(ut).gt.radii_fine(np,j).and.j.le.nfine)
           j = j + 1
          enddo
          partials_fineidx(ut) = j - 1
        enddo
      endif

!  Special case. Direct nadir viewing
!  ==================================

!  Compute everything and Exit.
!    (This is the same as the regular pseudo-spherical )

      if ( salpha_boa.eq.0.0d0 ) then

!  WHOLE LAYER and FINE divisions
!  ------------------------------

!  Start layer loop, working upwards

        do n = nlayers,1,-1

!  set main output.

          alpha_all(n-1)   = alpha_all(n)
          theta_all(n-1)   = theta_all(n)
          phi_all(n-1)     = phi_all(n)
          cosscat_up(n-1) = cosscat_up(n)
          lospaths(n) = radii(n-1)-radii(n)
          if ( do_fine ) then
            do j = 1, nfine
              alpha_fine(n,j) = 0.0d0
            enddo
          endif

!  Overhead sun

          if (stheta_boa.eq.0.0d0 ) then
            do k = n, 1, -1
              sunpaths(n,k) = radii(k-1)-radii(k)
            enddo
            if ( do_fine ) then
              do j = 1, nfine
                do k = n - 1, 1, -1
                  sunpaths_fine(n,k,j) = radii(k-1)-radii(k)
                enddo
                sunpaths_fine(n,n,j) = radii(n-1)-radii_fine(n,j)
              enddo
            endif
          endif

!  Non-overhead sun
!  Main output of solar paths
!  Solar path distances for fine output

          if (stheta_boa.gt.0.0d0 ) then
            sth0 = stheta_boa
            th0  = theta_all(n)
            do k = n, 1, -1
              sth1 = sth0*radii(k)/radii(k-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths(n,k) = dsin(ks1)*radii(k)/sth1
              sth0 = sth1
              th0  = th1
            enddo
            if ( do_fine ) then
              do j = 1, nfine
                sth0 = stheta_boa
                th0  = theta_all(n)
                sth1 = sth0*radii_fine(n,j)/radii(n-1)
                th1  = dasin(sth1)
                ks1  = th0-th1
                sunpaths_fine(n,n,j) = dsin(ks1)*radii_fine(n,j)/sth1
                sth0 = sth1
                th0  = th1
                do k = n-1, 1, -1
                  sth1 = sth0*radii(k)/radii(k-1)
                  th1  = dasin(sth1)
                  ks1  = th0-th1
                  sunpaths_fine(n,k,j) = dsin(ks1)*radii(k)/sth1
                  sth0 = sth1
                  th0  = th1
                enddo
              enddo
            endif
          endif

!  End main layer loop

        enddo

!  PARTIAL LAYERS
!  --------------

        if ( do_partials ) then
          do ut = 1, n_partials
            np = partials_idx(ut)

!  Los angle and paths

            alpha_p(ut) = alpha_all(nlayers)
!            lospaths_p_dn(ut) = radii(np-1)-radii_p(ut)
            lospaths_p_up(ut) = radii_p(ut)-radii(np)

!  Overhead sun

            if (stheta_boa.eq.0.0d0 ) then
              do k = 1, np-1
                sunpaths_p(ut,k) = radii(k-1)-radii(k)
              enddo
              sunpaths_p(ut,np) = radii(np-1)-radii_p(ut)
            endif

!  Non-overhead sun

            if (stheta_boa.gt.0.0d0 ) then
              sth0 = stheta_boa
              th0  = theta_all(np)
              sth1 = sth0*radii_p(ut)/radii(np-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths_p(ut,np) = dsin(ks1)*radii_p(ut)/sth1
              sth0 = sth1
              th0  = th1
              do k = np-1, 1, -1
                sth1 = sth0*radii(k)/radii(k-1)
                th1  = dasin(sth1)
                ks1  = th0-th1
                sunpaths_p(ut,k) = dsin(ks1)*radii(k)/sth1
                sth0 = sth1
                th0  = th1
              enddo
            endif

!  Finish partial layer stuff

          enddo
        endif

!  Return, as everything now done

        return

!  end regular pseudo-spherical clause, LOS is zero

      endif

!  Outgoing spehricity geometry
!  ============================

!  define Unit solar vector at BOA

      ex = - stheta_boa * cphi_boa
      ey = - stheta_boa * sphi_boa
      ez = - ctheta_boa

!  Sun paths, boa geometry, always directly illuminated

      if ( stheta_boa.eq.0.0d0 ) then
        do k = nlayers, 1, -1
          sunpaths(nlayers,k) = radii(k-1)-radii(k)
        enddo
      else
        sth0 = stheta_boa
        th0  = theta_all(nlayers)
        do k = nlayers, 1, -1
          sth1 = sth0*radii(k)/radii(k-1)
          th1  = dasin(sth1)
          ks1  = th0-th1
          sunpaths(nlayers,k) = dsin(ks1)*radii(k)/sth1
          sth0 = sth1
          th0  = th1
        enddo
      endif

!  Check single illumination
!      --Not required, now we have the tangent point treatment
!      if (stheta_boa.gt.0.0d0 ) then
!        xicum = dasin(radii(nlayers)*salpha_boa/radii(0))
!        xicum = alpha_all(nlayers)-xicum
!        px = - radii(0) * dsin(xicum)
!        py = 0.0d0
!        pz =   radii(0) * dcos(xicum)
!        b = ex*px + ey*py + ez*pz
!        ctheta = -b/radii(0)
!        if ( ctheta.le.0.0d0 ) then
!          write(*,*)'limit value = ',90.0d0-xicum/deg_to_rad
!        endif
!      endif

!  initialise los cumulative angle

      xicum  = 0.0d0

!  set TOA direct illumination flag

      direct_sun = .true.
      if ( do_fine ) then
        do j = 1, nfine
          direct_sunf(j) = .true.
        enddo
      endif

!  Start loop over positions (layer upper boundaries)

      do n = nlayers - 1, 0, -1

!  Next level up

        n1 = n + 1

!  Los angles at level boundaries

        salpha = radii(nlayers) * salpha_boa / radii(n)
        alpha_all(n)  = dasin(salpha)
        calpha = dcos(alpha_all(n))

!  Lospaths

        ksi = alpha_all(n1) - alpha_all(n)
        sksi = dsin(ksi)
        cksi = dcos(ksi)
        lospaths(n1) = sksi * radii(n1) / salpha
        xicum0 = xicum
        xicum  = xicum + ksi

!  Fine grid lospath output (angle and radius)
!    Locally save the earth-center angle ksif

        if ( do_fine ) then
          difa = (alpha_all(n1)-alpha_all(n))/dfine1
          do j = 1, nfine
            alpha_fine(n1,j) = alpha_all(n1) - difa * dble(j)
            saf = dsin(alpha_fine(n1,j))
            radii_fine(n1,j) = salpha_boa * radii(nlayers) / saf
            ksif(j) = alpha_all(n1) - alpha_fine(n1,j)
          enddo
        endif

!  Sun angles for the Direct Nadir case

        if (stheta_boa.eq.0.0d0 ) then
         theta_all(n) = xicum
         ctheta = dcos(theta_all(n))
         stheta = dsqrt(1.0d0-ctheta*ctheta)
         if ( do_fine ) then
           do j = 1, nfine
             thetaf(j)  = xicum0 + ksif(j)
             cthetaf(j) = dcos(thetaf(j))
             sthetaf(j) = dsqrt(1.0d0-ctheta*ctheta)
           enddo
         endif
        endif

!  Sun angles for the general case
!    Local save of angles, cosines, sines and  illumination flags

        if (stheta_boa.gt.0.0d0 ) then
         px = - radii(n) * dsin(xicum)
         py = 0.0d0
         pz =   radii(n) * dcos(xicum)
         b = ex*px + ey*py + ez*pz
         ctheta = -b/radii(n)
         direct_sun = (direct_sun.and.ctheta.ge.0.d0)
         stheta = dsqrt(1.0d0-ctheta*ctheta)
         theta_all(n) = dacos(ctheta)
         if ( do_fine ) then
           do j = 1, nfine
             xicumf  = xicum0 + ksif(j)
             px = - radii_fine(n1,j) * dsin(xicumf)
             py = 0.0d0
             pz =   radii_fine(n1,j) * dcos(xicumf)
             b  = ex*px + ey*py + ez*pz
             cthetaf(j) = -b/radii_fine(n1,j)
             direct_sunf(j) = (direct_sunf(j).and.cthetaf(j).ge.0.d0)
             sthetaf(j) = dsqrt(1.0d0-cthetaf(j)*cthetaf(j))
             thetaf(j)  = dacos(cthetaf(j))
           enddo
         endif
        endif

!  Unit vector f2(i) perpendicular to OP but in plane of path
!  projection of f2(i) on solar path gives the relative azimuth at P
!        f2x = dsin(xicum)
!        f2y = 0.0d0
!        f2z = dcos(xicum)
!        cphi = - (ex*f2x + ey*f2y + ez*f2z ) / stheta
!        cphi = - (ex*f2x + ey*f2y + ez*f2z ) / stheta
!        if ( cphi.gt.1.0d0)  cphi = 1.0d0
!        if ( cphi.lt.-1.0d0) cphi = -1.0d0
! ********************************************* Apparently not correct

!  Fix phi by using constancy of scatter angle

        cosscat_up(n) = cosscat_up(n+1)
        if (stheta_boa.eq.0.0d0 ) then
          phi_all(n)     = phi_all(n+1)
        else
         cphi = (cosscat_up(n)+calpha*ctheta)/stheta/salpha
         if ( cphi.gt.1.0d0) cphi = 1.0d0
         if ( cphi.lt.-1.0d0) cphi = -1.0d0
         phi_all(n)     = dacos(cphi)
         phi_all(n)     = dacos(cphi)
        endif

!  Sun paths, Direct sun at layer top
!  ==================================

!   Means that the SZA at layer top is < 90.
!    ===> SZA < 90 for all fine points in layer beneath

        if ( direct_sun ) then

!  Work up from level n to TOA
!    Layer top calculation gets left out at TOA

         if ( n .gt. 0 ) then
          sth0 = stheta
          th0  = theta_all(n)
          do k = n, 1, -1
           sth1 = sth0*radii(k)/radii(k-1)
           th1  = dasin(sth1)
           ks1  = th0-th1
           sunpaths(n,k) = dsin(ks1)*radii(k)/sth1
           sth0 = sth1
           th0  = th1
          enddo
         endif

! DBG         write(*,*)'regular',n1,(sunpaths(n1,k),k=n1,1,-1)

!  Fine grid calculation is always required
!  ----------------------------------------

!   Start at the grid point on LOS path and work upwards,
!     first across partial layer to upper boundary, then continue
!     upwards to TOA on a whole layer basis. Sine rule.

         if ( do_fine ) then
          do j = 1, nfine
           sth0 = sthetaf(j)
           th0  = thetaf(j)
           sth1 = sth0*radii_fine(n1,j)/radii(n)
           th1  = dasin(sth1)
           ks1  = th0-th1
           sunpaths_fine(n1,n1,j) = dsin(ks1)*radii_fine(n1,j)/sth1
           sth0 = sth1
           th0  = th1
           do k = n, 1, -1
            sth1 = sth0*radii(k)/radii(k-1)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k)/sth1
            sth0 = sth1
            th0  = th1
           enddo
!  DBG           write(*,*)'regular',n1,j,(sunpaths_fine(n1,k,j),k=n1,1,-1)
          enddo

!  DBG         write(*,*)'regular',n,(sunpaths(n,k),k=n,1,-1)

         endif

!  Complete direct sun computations

        endif

!  Sun paths, Not direct sun , with tangent point
!  ==============================================

!  Although layer top has a tangent point, not all of the fine-grid
!   points will have a tangent point.

        if (.not.direct_sun ) then

!  First do the layer-top calculation.
!  -----------------------------------

!  TANGR = tangent point radius.

         tangr = stheta*radii(n)

!  ntraverse(n) is the number of layers traversed by ray.

         krad = nlayers
         do while (tangr.gt.radii(krad))
          krad = krad - 1
         enddo
         ntraverse(n) = krad + 1

!  Start at the TOA angles

         sth0 = tangr/radii(0)
         th0 = dasin(sth0)

!  Work downwards from TOA (sine rule) to level immediately above
!  the tangent layer. Don't forget to double the path length for
!  any layers which are traversed twice.

         do k = 1, krad
          sth1 = radii(k-1)*sth0/radii(k)
          th1 = dasin(sth1)
          ks1 = th1-th0
          fac = 1.0d0
          if ( k.gt.n) fac = 2.0d0
          sunpaths(n,k) = fac*dsin(ks1)*radii(k-1)/sth1
          sth0 = sth1
          th0  = th1
         enddo

!  Tangent layer path length. Twice again. The following check is good.
!  check       write(*,*)tangr/dtan(th1),radii(krad)*dcos(th1)

         sunpaths(n,krad+1)=2.0d0*radii(krad)*dcos(th1)

! DBG         write(*,*)'--tangent',n1,(sunpaths(n1,k),k=ntraverse(n1),1,-1)

!  Fine layer development (slightly different)
!  ---------------------

         if ( do_fine ) then
          do j = 1, nfine

!  If there is no tangent point, repeat calculation above.

           if ( direct_sunf(j) ) then
            sth0 = sthetaf(j)
            th0  = thetaf(j)
            sth1 = sth0*radii_fine(n1,j)/radii(n)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_fine(n1,n1,j) = dsin(ks1)*radii_fine(n1,j)/sth1
            sth0 = sth1
            th0  = th1
            do k = n, 1, -1
             sth1 = sth0*radii(k)/radii(k-1)
             th1  = dasin(sth1)
             ks1  = th0-th1
             sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k)/sth1
             sth0 = sth1
             th0  = th1
            enddo

! DBG           write(*,*)'regular',n1,j,
!     &       (sunpaths_fine(n1,k,j),k=ntraverse_fine(n1,j),1,-1)

!  Fine grid Calculation with tangent point
!  ----------------------------------------

           else

!  Local tangent radius and number of layers traversed

            tangr = sthetaf(j)*radii_fine(n1,j)
            krad = nlayers
            do while (tangr.gt.radii(krad))
             krad = krad - 1
            enddo
            ntraverse_fine(n1,j) = krad + 1

!  Start again at TOA

            sth0 = tangr/radii(0)
            th0  = dasin(sth0)

!  Work down to the level n

            do k = 1, n
             sth1 = radii(k-1)*sth0/radii(k)
             th1 = dasin(sth1)
             ks1 = th1-th0
             sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k-1)/sth1
             sth0 = sth1
             th0  = th1
            enddo

!  In layer below level n, Work down to level of fine grid radius
!     (single contribution)

            sth1 = radii(n)*sth0/radii_fine(n1,j)
            th1 = dasin(sth1)
            ks1 = th1-th0
            path = dsin(ks1)*radii(n)/sth1
            sth0 = sth1
            th0  = th1

!  In layer below level n, Complete the path down to the tangent point
!    (double contribution). Finish.

            if ( krad.eq.n ) then

              path = path + 2.0d0*radii_fine(n1,j)*dcos(th1)
              sunpaths_fine(n1,krad+1,j) = path

!  Check    write(*,*)'1',tangr/dtan(th1),radii_fine(n1,j)*dcos(th1)

!  In layer below level n, Need to go down further.
!    (from now on, we need double contributions)
!     --- Continue the path to the bottom of this layer
!     --- Continue down by whole-layer steps until reach level above Tangent
!     --- Complete the path down to the tangent point

            else
              sth1 = radii_fine(n1,j)*sth0/radii(n1)
              th1 = dasin(sth1)
              ks1 = th1-th0
              path = path + 2.0d0 * dsin(ks1)*radii_fine(n1,j)/sth1
              sunpaths_fine(n1,n1,j) = path
              sth0 = sth1
              th0  = th1
              do k = n1 + 1, krad
               sth1 = radii(k-1)*sth0/radii(k)
               th1 = dasin(sth1)
               ks1 = th1-th0
               sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k-1)/sth1
               sth0 = sth1
               th0  = th1
              enddo
              sunpaths_fine(n1,krad+1,j)=2.0d0*radii(krad)*dcos(th1)
!  Check 2    write(*,*)'2',tangr/dtan(th1),radii(krad)*dcos(th1)
            endif

! DBG           write(*,*)'tangent',n1,j,
!     &       (sunpaths_fine(n1,k,j),k=ntraverse_fine(n1,j),1,-1)

!  Complete tangent point clause

           endif

!  Complete fine-grid loop

          enddo

!  DBG        write(*,*)'tangent',n,(sunpaths(n,k),k=ntraverse(n),1,-1)

!  Complete tangent point calculation

         endif

        endif

!  End layer loop

      enddo

!  PARTIALS Section
!  ----------------

!  Partial angles and lospaths

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)

!  establish the LOS angle, and lospaths

          th0 = alpha_all(np)
          sth0 = dsin(th0)
          salpha = radii(np)*sth0 / radii_p(ut)
          alpha_p(ut) = dasin(salpha)
          calpha = dsqrt(1.0d0-salpha*salpha)

          ksi = alpha_all(np) - alpha_p(ut)
          lospaths_p_up(ut) = dsin(ksi)*radii(np)/salpha

!          ksi = alpha_p(ut) - alpha_all(np-1)
!          lospaths_p_dn(ut) = dsin(ksi)*radii(np-1)/salpha

!  Check debug
!          lospaths_p_dn(ut) = lospaths(np)-lospaths_p_up(ut)
!          write(*,*)lospaths_p_up(ut),lospaths_p_dn(ut),
!     &              radii_p(ut)-radii(np)
!          pause

!  Establish the solar angle, both naidr =sun-at-BOA and generally

          xicum = alpha_all(nlayers) - alpha_p(ut)
          direct_sun = .true.
          if (stheta_boa.eq.0.0d0 ) then
            theta_p = xicum
            ctheta = dcos(theta_p)
            stheta = dsqrt(1.0d0-ctheta*ctheta)
          else
            px = - radii_p(ut) * dsin(xicum)
            py = 0.0d0
            pz =   radii_p(ut) * dcos(xicum)
            b = ex*px + ey*py + ez*pz
            ctheta = -b/radii_p(ut)
            direct_sun = (direct_sun.and.ctheta.ge.0.d0)
            stheta = dsqrt(1.0d0-ctheta*ctheta)
            theta_p = dacos(ctheta)
          endif

!         write(*,*)direct_sun,theta_all(Np-1),theta_p,theta_all(np)

!  Sun paths, Direct sun at layer top
!   First calculate sunpath for layer containing partial point
!   Then work upwards to TOA using the sine rule.

          if ( direct_sun ) then
            sth0 = stheta
            th0  = theta_p
            sth1 = sth0*radii_p(ut)/radii(np-1)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_p(ut,np) = dsin(ks1)*radii_p(ut)/sth1
            sth0 = sth1
            th0  = th1
            do k = np-1, 1, -1
              sth1 = sth0*radii(k)/radii(k-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths_p(ut,k) = dsin(ks1)*radii(k)/sth1
              sth0 = sth1
              th0  = th1
            enddo
          endif

!  Debug
! 14       format(a11,1x,10f10.6)
!          j = partials_fineidx(ut)
!          write(*,*)ut,ntraverse_p(ut)
!          write(*,*)ut,j,np,lospaths_p_up(ut)/lospaths(np)
!           write(*,14)'regular feb',(sunpaths(np,k),k=np,1,-1)
!           do j = 1, nfine
!             write(*,14)'regular feb',(sunpaths_fine(np,k,j),k=np,1,-1)
!             if (j.eq.partials_fineidx(ut))
!     &        write(*,14)'regular ut ',(sunpaths_p(ut,k),k=np,1,-1)
!           enddo
!           write(*,14)'regular feb',0.0d0,(sunpaths(np-1,k),k=np-1,1,-1)
!
!  Sun paths, Not direct sun , with tangent point. Code similar to above.
!  TANGR = tangent point radius for this ray
!   ntraverse_p(ut) = Number of layers traversed by ray.

!  Work downwards from TOA (sine rule) to level immediately above
!  the tangent layer. Don't forget to double the path length for
!  any layers which are traversed twice.

!  Tangent layer path length. Twice again. The following check is good.
!  check       write(*,*)tangr/dtan(th1),radii(krad)*dcos(th1)

          if (.not.direct_sun ) then
            tangr = stheta*radii(n)
            krad = nlayers
            do while (tangr.gt.radii(krad))
              krad = krad - 1
            enddo
            ntraverse_p(ut) = krad + 1
            sth0 = tangr/radii(0)
            th0 = dasin(sth0)
            do k = 1, krad
              sth1 = radii(k-1)*sth0/radii(k)
              th1 = dasin(sth1)
              ks1 = th1-th0
              fac = 1.0d0
              if ( k.gt.np) fac = 2.0d0
              sunpaths_p(ut,k) = fac*dsin(ks1)*radii(k-1)/sth1
              sth0 = sth1
              th0  = th1
            enddo
            sunpaths_p(ut,krad+1)=2.0d0*radii_p(ut)*dcos(th1)
          endif

!  Finish partials loop

        enddo
      endif

!  Finish

      return
END SUBROUTINE outgoing_sphergeom_fine_up

!

SUBROUTINE outgoing_sphergeom_fine_dn                              &
       ( maxlayers, maxfine, maxpartials,                          & ! Input
         do_fine, do_partials, nlayers, nfine,                     & ! Input
         n_partials, partials_idx,                                 & ! Input
         heights, heights_p, eradius,                              & ! Input
         alpha_boa, theta_boa, phi_boa,                            & ! Input
         sunpaths,      radii,      ntraverse,      alpha_all,     & ! Output
         sunpaths_fine, ntraverse_fine, alpha_fine,                & ! Output
         sunpaths_p,    ntraverse_p,    alpha_p,                   & ! Output
         partials_fineidx, cosscat_dn, fail, message )               ! Output

!  Completely stand-alone geometry routine for the outgoing correction
!    This is applicable to the Downwelling path geoemtry

!    starting inputs are the BOA values of SZA, VZA and PHI
!    need also the height grids, earth radius and control

!  This routine has the fine gridding treatment
!  Version 2.3. September 2007, Partial layer geometries added

      implicit none

!  inputs

      integer       :: maxlayers, maxfine, maxpartials
      integer       :: nlayers, nfine
      logical       :: do_fine, do_partials
      integer       :: n_partials
      integer       :: partials_idx    (maxpartials)
      real(fpk)     :: eradius, heights (0:maxlayers)
      real(fpk)     :: heights_p(maxpartials)
      real(fpk)     :: alpha_boa, theta_boa, phi_boa

!  main outputs (geometry)

      integer       :: ntraverse   (0:maxlayers)
      real(fpk)     :: sunpaths   (0:maxlayers,maxlayers)
      real(fpk)     :: radii      (0:maxlayers)
      real(fpk)     :: alpha_all  (0:maxlayers)

!  Fine level output (geometry)

      integer       :: ntraverse_fine(maxlayers,maxfine)
      real(fpk)     :: sunpaths_fine (maxlayers,maxlayers,maxfine)
      real(fpk)     :: alpha_fine    (maxlayers,maxfine)

!  Partial layer output (geometry)

      integer       :: ntraverse_p(maxpartials)
      real(fpk)     :: sunpaths_p (maxpartials,maxlayers)
      real(fpk)     :: alpha_p    (maxpartials)

!  scattering angles

      real(fpk)     :: cosscat_dn (0:maxlayers)

!   Local arrays

      real(fpk)     :: radii_fine    (maxlayers,maxfine)
      real(fpk)     :: radii_p    (maxpartials)
      real(fpk)     :: lospaths(maxlayers)

      real(fpk)     :: lospaths_p_dn(maxpartials)
      integer       :: partials_fineidx(maxpartials)

      real(fpk)     :: theta_all  (0:maxlayers)
      real(fpk)     :: phi_all    (0:maxlayers)

!  Status output

      logical       :: fail
      character*(*) :: message

!  Local

      logical       :: direct_sun
      integer       :: n, j, k, krad, n1, ut, np
      real(fpk)     :: deg_to_rad, ex, ey, ez, px, py, pz
      real(fpk)     :: salpha_boa, calpha_boa, sphi_boa
      real(fpk)     :: stheta_boa, ctheta_boa, cphi_boa
      real(fpk)     :: ksi, cksi, sksi, xicum, tangr, fac
      real(fpk)     :: ctheta, stheta, calpha, salpha, cphi
      real(fpk)     :: b, sth0, th0, ks1, sth1, th1, theta_p

!  Local arrays associated with fine grid output

      integer, parameter :: maxlocalfine = 20
      logical            :: direct_sunf(maxlocalfine)
      real(fpk)          :: difz, dfine1, saf, xicum0, path
      real(fpk)          :: thetaf(maxlocalfine), xicumf, difa
      real(fpk)          :: cthetaf(maxlocalfine)
      real(fpk)          :: sthetaf(maxlocalfine)
      real(fpk)          :: ksif(maxlocalfine)

!  Initialise output

      fail = .false.
      message = ' '

!  check range of inputs

      if ( alpha_boa.ge.90.0d0.or.alpha_boa.lt.0.0d0 ) then
        message = 'boa LOS angle outside range [0,90]); Check it!'
        fail    = .true.
        return
      endif
      if ( phi_boa.lt.0.0d0 )   phi_boa = - phi_boa
      if ( phi_boa.gt.180.0d0 ) phi_boa = 360.0d0 - phi_boa
      if ( theta_boa.ge.90.0d0.or.theta_boa.lt.0.0d0 ) then
        message = 'boa SZA angle outside range [0,90]); Check it!'
        fail    = .true.
        return
      endif
      if ( do_fine ) then
       if ( nfine.gt.maxlocalfine ) then
         message = 'local finelayer dimensioning insufficient'
         fail    = .true.
         return
       endif
      endif

!  zero the sun paths
!  Initialize number of layers traversed  (nominal conditions)

      do n = 0, nlayers
        ntraverse(n) = n
        do k = 1, nlayers
         sunpaths(n,k) = 0.0d0
        enddo
      enddo

!  Zero the partial paths, initialize the traverse number (nominal)

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)
          ntraverse_p(ut) = np
          do k = 1, np
            sunpaths_p(ut,k) = 0.0d0
          enddo
        enddo
      endif

!  Zero the fine data paths

      if ( do_fine ) then
       dfine1 = dble(nfine) + 1
       do n = 1, nlayers
        do j = 1, nfine
         ntraverse_fine(n,j) = n
         do k = 1, nlayers
          sunpaths_fine(n,k,j) = 0.0d0
         enddo
        enddo
       enddo
      endif

!  start at BOA

      deg_to_rad = dacos(-1.0d0) / 180.0d0
      alpha_all(nlayers) = alpha_boa * deg_to_rad
      theta_all(nlayers) = theta_boa * deg_to_rad
      phi_all(nlayers)   = phi_boa   * deg_to_rad

!  Cosine of scattering angle at boa

      salpha_boa = dsin(alpha_all(nlayers))
      calpha_boa = dcos(alpha_all(nlayers))
      stheta_boa = dsin(theta_all(nlayers))
      ctheta_boa = dcos(theta_all(nlayers))
      cphi_boa   = dcos(phi_all(nlayers))
      sphi_boa   = dsin(phi_all(nlayers))

      cosscat_dn (nlayers) = + calpha_boa * ctheta_boa + &
                               salpha_boa * stheta_boa * cphi_boa

!  Radii
!  -----

!  layer levels

      do n = 0, nlayers
        radii(n) = eradius + heights(n)
      enddo

!  Fine levels

      if ( do_fine ) then
        do n = 1, nlayers
          difz = (radii(n-1)-radii(n))/dfine1
          do j = 1, nfine
            radii_fine(n,j) = radii(n) + difz * dble(j)
          enddo
        enddo
      endif

!  Partial levels
!   Find the fine level just below partial point = (0,1,2,...nfine)

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)
          radii_p(ut) = eradius + heights_p(ut)
          j = 1
          do while (radii_p(ut).gt.radii_fine(np,j).and.j.le.nfine)
           j = j + 1
          enddo
          partials_fineidx(ut) = j - 1
        enddo
      endif

!  Special case. Direct nadir viewing
!  ==================================

!  Compute everything and Exit.
!    (This is the same as the regular pseudo-spherical )

      if ( salpha_boa.eq.0.0d0 ) then

!  WHOLE LAYER and FINE divisions
!  ------------------------------

!  Start layer loop, working upwards

        do n = nlayers,1,-1

!  set main output.

          alpha_all(n-1)   = alpha_all(n)
          theta_all(n-1)   = theta_all(n)
          phi_all(n-1)     = phi_all(n)
          cosscat_dn(n-1) = cosscat_dn(n)
          lospaths(n) = radii(n-1)-radii(n)
          if ( do_fine ) then
            do j = 1, nfine
              alpha_fine(n,j) = 0.0d0
            enddo
          endif

!  Overhead sun

          if (stheta_boa.eq.0.0d0 ) then
            do k = n, 1, -1
              sunpaths(n,k) = radii(k-1)-radii(k)
            enddo
            if ( do_fine ) then
              do j = 1, nfine
                do k = n - 1, 1, -1
                  sunpaths_fine(n,k,j) = radii(k-1)-radii(k)
                enddo
                sunpaths_fine(n,n,j) = radii(n-1)-radii_fine(n,j)
              enddo
            endif
          endif

!  Non-overhead sun
!  Main output of solar paths
!  Solar path distances for fine output

          if (stheta_boa.gt.0.0d0 ) then
            sth0 = stheta_boa
            th0  = theta_all(n)
            do k = n, 1, -1
              sth1 = sth0*radii(k)/radii(k-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths(n,k) = dsin(ks1)*radii(k)/sth1
              sth0 = sth1
              th0  = th1
            enddo
            if ( do_fine ) then
              do j = 1, nfine
                sth0 = stheta_boa
                th0  = theta_all(n)
                sth1 = sth0*radii_fine(n,j)/radii(n-1)
                th1  = dasin(sth1)
                ks1  = th0-th1
                sunpaths_fine(n,n,j) = dsin(ks1)*radii_fine(n,j)/sth1
                sth0 = sth1
                th0  = th1
                do k = n-1, 1, -1
                  sth1 = sth0*radii(k)/radii(k-1)
                  th1  = dasin(sth1)
                  ks1  = th0-th1
                  sunpaths_fine(n,k,j) = dsin(ks1)*radii(k)/sth1
                  sth0 = sth1
                  th0  = th1
                enddo
              enddo
            endif
          endif

!  End main layer loop

        enddo

!  PARTIAL LAYERS
!  --------------

        if ( do_partials ) then
          do ut = 1, n_partials
            np = partials_idx(ut)

!  Los angle and paths

            alpha_p(ut) = alpha_all(nlayers)
            lospaths_p_dn(ut) = radii(np-1)-radii_p(ut)
!            lospaths_p_up(ut) = radii_p(ut)-radii(np)

!  Overhead sun

            if (stheta_boa.eq.0.0d0 ) then
              do k = 1, np-1
                sunpaths_p(ut,k) = radii(k-1)-radii(k)
              enddo
              sunpaths_p(ut,np) = radii(np-1)-radii_p(ut)
            endif

!  Non-overhead sun

            if (stheta_boa.gt.0.0d0 ) then
              sth0 = stheta_boa
              th0  = theta_all(np)
              sth1 = sth0*radii_p(ut)/radii(np-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths_p(ut,np) = dsin(ks1)*radii_p(ut)/sth1
              sth0 = sth1
              th0  = th1
              do k = np-1, 1, -1
                sth1 = sth0*radii(k)/radii(k-1)
                th1  = dasin(sth1)
                ks1  = th0-th1
                sunpaths_p(ut,k) = dsin(ks1)*radii(k)/sth1
                sth0 = sth1
                th0  = th1
              enddo
            endif

!  Finish partial layer stuff

          enddo
        endif

!  Return, as everything now done

        return

!  end regular pseudo-spherical clause, LOS is zero

      endif

!  Outgoing sphericity geometry
!  ============================

!  define Unit solar vector at BOA
!   Note change of sign for ex here.

      ex = + stheta_boa * cphi_boa
      ey = - stheta_boa * sphi_boa
      ez = - ctheta_boa

!  Sun paths, boa geometry, always directly illuminated

      if ( stheta_boa.eq.0.0d0 ) then
        do k = nlayers, 1, -1
          sunpaths(nlayers,k) = radii(k-1)-radii(k)
        enddo
      else
        sth0 = stheta_boa
        th0  = theta_all(nlayers)
        do k = nlayers, 1, -1
          sth1 = sth0*radii(k)/radii(k-1)
          th1  = dasin(sth1)
          ks1  = th0-th1
          sunpaths(nlayers,k) = dsin(ks1)*radii(k)/sth1
          sth0 = sth1
          th0  = th1
        enddo
      endif

!  Check single illumination
!      --Not required, now we have the tangent point treatment
!      if (stheta_boa.gt.0.0d0 ) then
!        xicum = dasin(radii(nlayers)*salpha_boa/radii(0))
!        xicum = alpha_all(nlayers)-xicum
!        px = - radii(0) * dsin(xicum)
!        py = 0.0d0
!        pz =   radii(0) * dcos(xicum)
!        b = ex*px + ey*py + ez*pz
!        ctheta = -b/radii(0)
!        if ( ctheta.le.0.0d0 ) then
!          write(*,*)'limit value = ',90.0d0-xicum/deg_to_rad
!        endif
!      endif

!  initialise los cumulative angle

      xicum  = 0.0d0

!  set TOA direct illumination flag

      direct_sun = .true.
      if ( do_fine ) then
        do j = 1, nfine
          direct_sunf(j) = .true.
        enddo
      endif

!  Start loop over positions (layer upper boundaries)

      do n = nlayers - 1, 0, -1

!  Next level up

        n1 = n + 1

!  Los angles at level boundaries

        salpha = radii(nlayers) * salpha_boa / radii(n)
        alpha_all(n)  = dasin(salpha)
        calpha = dcos(alpha_all(n))

!  Lospaths

        ksi = alpha_all(n1) - alpha_all(n)
        sksi = dsin(ksi)
        cksi = dcos(ksi)
        lospaths(n1) = sksi * radii(n1) / salpha
        xicum0 = xicum
        xicum  = xicum + ksi

!  Fine grid lospath output (angle and radius)
!    Locally save the earth-center angle ksif

        if ( do_fine ) then
          difa = (alpha_all(n1)-alpha_all(n))/dfine1
          do j = 1, nfine
            alpha_fine(n1,j) = alpha_all(n1) - difa * dble(j)
            saf = dsin(alpha_fine(n1,j))
            radii_fine(n1,j) = salpha_boa * radii(nlayers) / saf
            ksif(j) = alpha_all(n1) - alpha_fine(n1,j)
          enddo
        endif

!  Sun angles for the Direct Nadir case

        if (stheta_boa.eq.0.0d0 ) then
         theta_all(n) = xicum
         ctheta = dcos(theta_all(n))
         stheta = dsqrt(1.0d0-ctheta*ctheta)
         if ( do_fine ) then
           do j = 1, nfine
             thetaf(j)  = xicum0 + ksif(j)
             cthetaf(j) = dcos(thetaf(j))
             sthetaf(j) = dsqrt(1.0d0-ctheta*ctheta)
           enddo
         endif
        endif

!  Sun angles for the general case
!    Local save of angles, cosines, sines and  illumination flags

        if (stheta_boa.gt.0.0d0 ) then
         px = - radii(n) * dsin(xicum)
         py = 0.0d0
         pz =   radii(n) * dcos(xicum)
         b = ex*px + ey*py + ez*pz
         ctheta = -b/radii(n)
         direct_sun = (direct_sun.and.ctheta.ge.0.d0)
         stheta = dsqrt(1.0d0-ctheta*ctheta)
         theta_all(n) = dacos(ctheta)
         if ( do_fine ) then
           do j = 1, nfine
             xicumf  = xicum0 + ksif(j)
             px = - radii_fine(n1,j) * dsin(xicumf)
             py = 0.0d0
             pz =   radii_fine(n1,j) * dcos(xicumf)
             b  = ex*px + ey*py + ez*pz
             cthetaf(j) = -b/radii_fine(n1,j)
             direct_sunf(j) = (direct_sunf(j).and.cthetaf(j).ge.0.d0)
             sthetaf(j) = dsqrt(1.0d0-cthetaf(j)*cthetaf(j))
             thetaf(j)  = dacos(cthetaf(j))
           enddo
         endif
        endif

!  Unit vector f2(i) perpendicular to OP but in plane of path
!  projection of f2(i) on solar path gives the relative azimuth at P
!        f2x = dsin(xicum)
!        f2y = 0.0d0
!        f2z = dcos(xicum)
!        cphi = - (ex*f2x + ey*f2y + ez*f2z ) / stheta
!        cphi = - (ex*f2x + ey*f2y + ez*f2z ) / stheta
!        if ( cphi.gt.1.0d0)  cphi = 1.0d0
!        if ( cphi.lt.-1.0d0) cphi = -1.0d0
! ********************************************* Apparently not correct

!  Fix phi by using constancy of scatter angle

        cosscat_dn(n) = cosscat_dn(n+1)
        if (stheta_boa.eq.0.0d0 ) then
          phi_all(n)     = phi_all(n+1)
        else
         cphi = (cosscat_dn(n)-calpha*ctheta)/stheta/salpha
         if ( cphi.gt.1.0d0) cphi = 1.0d0
         if ( cphi.lt.-1.0d0) cphi = -1.0d0
         phi_all(n)     = dacos(cphi)
         phi_all(n)     = dacos(cphi)
        endif

!  Sun paths, Direct sun at layer top
!  ==================================

!   Means that the SZA at layer top is < 90.
!    ===> SZA < 90 for all fine points in layer beneath

        if ( direct_sun ) then

!  Work up from level n to TOA
!    Layer top calculation gets left out at TOA

         if ( n .gt. 0 ) then
          sth0 = stheta
          th0  = theta_all(n)
          do k = n, 1, -1
           sth1 = sth0*radii(k)/radii(k-1)
           th1  = dasin(sth1)
           ks1  = th0-th1
           sunpaths(n,k) = dsin(ks1)*radii(k)/sth1
           sth0 = sth1
           th0  = th1
          enddo
         endif

! DBG         write(*,*)'regular',n1,(sunpaths(n1,k),k=n1,1,-1)

!  Fine grid calculation is always required
!  ----------------------------------------

!   Start at the grid point on LOS path and work upwards,
!     first across partial layer to upper boundary, then continue
!     upwards to TOA on a whole layer basis. Sine rule.

         if ( do_fine ) then
          do j = 1, nfine
           sth0 = sthetaf(j)
           th0  = thetaf(j)
           sth1 = sth0*radii_fine(n1,j)/radii(n)
           th1  = dasin(sth1)
           ks1  = th0-th1
           sunpaths_fine(n1,n1,j) = dsin(ks1)*radii_fine(n1,j)/sth1
           sth0 = sth1
           th0  = th1
           do k = n, 1, -1
            sth1 = sth0*radii(k)/radii(k-1)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k)/sth1
            sth0 = sth1
            th0  = th1
           enddo
!  DBG           write(*,*)'regular',n1,j,(sunpaths_fine(n1,k,j),k=n1,1,-1)
          enddo

!  DBG         write(*,*)'regular',n,(sunpaths(n,k),k=n,1,-1)

         endif

!  Complete direct sun computations

        endif

!  Sun paths, Not direct sun , with tangent point
!  ==============================================

!  Although layer top has a tangent point, not all of the fine-grid
!   points will have a tangent point.

        if (.not.direct_sun ) then

!  First do the layer-top calculation.
!  -----------------------------------

!  TANGR = tangent point radius.

         tangr = stheta*radii(n)

!  ntraverse(n) is the number of layers traversed by ray.

         krad = nlayers
         do while (tangr.gt.radii(krad))
          krad = krad - 1
         enddo
         ntraverse(n) = krad + 1

!  Start at the TOA angles

         sth0 = tangr/radii(0)
         th0 = dasin(sth0)

!  Work downwards from TOA (sine rule) to level immediately above
!  the tangent layer. Don't forget to double the path length for
!  any layers which are traversed twice.

         do k = 1, krad
          sth1 = radii(k-1)*sth0/radii(k)
          th1 = dasin(sth1)
          ks1 = th1-th0
          fac = 1.0d0
          if ( k.gt.n) fac = 2.0d0
          sunpaths(n,k) = fac*dsin(ks1)*radii(k-1)/sth1
          sth0 = sth1
          th0  = th1
         enddo

!  Tangent layer path length. Twice again. The following check is good.
!  check       write(*,*)tangr/dtan(th1),radii(krad)*dcos(th1)

         sunpaths(n,krad+1)=2.0d0*radii(krad)*dcos(th1)

! DBG         write(*,*)'--tangent',n1,(sunpaths(n1,k),k=ntraverse(n1),1,-1)

!  Fine layer development (slightly different)
!  ---------------------

         if ( do_fine ) then
          do j = 1, nfine

!  If there is no tangent point, repeat calculation above.

           if ( direct_sunf(j) ) then
            sth0 = sthetaf(j)
            th0  = thetaf(j)
            sth1 = sth0*radii_fine(n1,j)/radii(n)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_fine(n1,n1,j) = dsin(ks1)*radii_fine(n1,j)/sth1
            sth0 = sth1
            th0  = th1
            do k = n, 1, -1
             sth1 = sth0*radii(k)/radii(k-1)
             th1  = dasin(sth1)
             ks1  = th0-th1
             sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k)/sth1
             sth0 = sth1
             th0  = th1
            enddo

! DBG           write(*,*)'regular',n1,j,
!     &       (sunpaths_fine(n1,k,j),k=ntraverse_fine(n1,j),1,-1)

!  Fine grid Calculation with tangent point
!  ----------------------------------------

           else

!  Local tangent radius and number of layers traversed

            tangr = sthetaf(j)*radii_fine(n1,j)
            krad = nlayers
            do while (tangr.gt.radii(krad))
             krad = krad - 1
            enddo
            ntraverse_fine(n1,j) = krad + 1

!  Start again at TOA

            sth0 = tangr/radii(0)
            th0  = dasin(sth0)

!  Work down to the level n

            do k = 1, n
             sth1 = radii(k-1)*sth0/radii(k)
             th1 = dasin(sth1)
             ks1 = th1-th0
             sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k-1)/sth1
             sth0 = sth1
             th0  = th1
            enddo

!  In layer below level n, Work down to level of fine grid radius
!     (single contribution)

            sth1 = radii(n)*sth0/radii_fine(n1,j)
            th1 = dasin(sth1)
            ks1 = th1-th0
            path = dsin(ks1)*radii(n)/sth1
            sth0 = sth1
            th0  = th1

!  In layer below level n, Complete the path down to the tangent point
!    (double contribution). Finish.

            if ( krad.eq.n ) then

              path = path + 2.0d0*radii_fine(n1,j)*dcos(th1)
              sunpaths_fine(n1,krad+1,j) = path

!  Check    write(*,*)'1',tangr/dtan(th1),radii_fine(n1,j)*dcos(th1)

!  In layer below level n, Need to go down further.
!    (from now on, we need double contributions)
!     --- Continue the path to the bottom of this layer
!     --- Continue down by whole-layer steps until reach level above Tangent
!     --- Complete the path down to the tangent point

            else
              sth1 = radii_fine(n1,j)*sth0/radii(n1)
              th1 = dasin(sth1)
              ks1 = th1-th0
              path = path + 2.0d0 * dsin(ks1)*radii_fine(n1,j)/sth1
              sunpaths_fine(n1,n1,j) = path
              sth0 = sth1
              th0  = th1
              do k = n1 + 1, krad
               sth1 = radii(k-1)*sth0/radii(k)
               th1 = dasin(sth1)
               ks1 = th1-th0
               sunpaths_fine(n1,k,j) = dsin(ks1)*radii(k-1)/sth1
               sth0 = sth1
               th0  = th1
              enddo
              sunpaths_fine(n1,krad+1,j)=2.0d0*radii(krad)*dcos(th1)
!  Check 2    write(*,*)'2',tangr/dtan(th1),radii(krad)*dcos(th1)
            endif

! DBG           write(*,*)'tangent',n1,j,
!     &       (sunpaths_fine(n1,k,j),k=ntraverse_fine(n1,j),1,-1)

!  Complete tangent point clause

           endif

!  Complete fine-grid loop

          enddo

!  DBG        write(*,*)'tangent',n,(sunpaths(n,k),k=ntraverse(n),1,-1)

!  Complete tangent point calculation

         endif

        endif

!  End layer loop

      enddo

!  PARTIALS Section
!  ----------------

!  Partial angles and lospaths

      if ( do_partials ) then
        do ut = 1, n_partials
          np = partials_idx(ut)

!  establish the LOS angle, and lospaths

          th0 = alpha_all(np)
          sth0 = dsin(th0)
          salpha = radii(np)*sth0 / radii_p(ut)
          alpha_p(ut) = dasin(salpha)
          calpha = dsqrt(1.0d0-salpha*salpha)

!          ksi = alpha_all(np) - alpha_p(ut)
!          lospaths_p_up(ut) = dsin(ksi)*radii(np)/salpha

          ksi = alpha_p(ut) - alpha_all(np-1)
          lospaths_p_dn(ut) = dsin(ksi)*radii(np-1)/salpha

!  Establish the solar angle, both naidr =sun-at-BOA and generally

          xicum = alpha_all(nlayers) - alpha_p(ut)
          direct_sun = .true.
          if (stheta_boa.eq.0.0d0 ) then
            theta_p = xicum
            ctheta = dcos(theta_p)
            stheta = dsqrt(1.0d0-ctheta*ctheta)
          else
            px = - radii_p(ut) * dsin(xicum)
            py = 0.0d0
            pz =   radii_p(ut) * dcos(xicum)
            b = ex*px + ey*py + ez*pz
            ctheta = -b/radii_p(ut)
            direct_sun = (direct_sun.and.ctheta.ge.0.d0)
            stheta = dsqrt(1.0d0-ctheta*ctheta)
            theta_p = dacos(ctheta)
          endif

!         write(*,*)direct_sun,theta_all(Np-1),theta_p,theta_all(np)

!  Sun paths, Direct sun at layer top
!   First calculate sunpath for layer containing partial point
!   Then work upwards to TOA using the sine rule.

          if ( direct_sun ) then
            sth0 = stheta
            th0  = theta_p
            sth1 = sth0*radii_p(ut)/radii(np-1)
            th1  = dasin(sth1)
            ks1  = th0-th1
            sunpaths_p(ut,np) = dsin(ks1)*radii_p(ut)/sth1
            sth0 = sth1
            th0  = th1
            do k = np-1, 1, -1
              sth1 = sth0*radii(k)/radii(k-1)
              th1  = dasin(sth1)
              ks1  = th0-th1
              sunpaths_p(ut,k) = dsin(ks1)*radii(k)/sth1
              sth0 = sth1
              th0  = th1
            enddo
          endif

!  Debug
! 14       format(a11,1x,10f10.6)
!          j = partials_fineidx(ut)
!          write(*,*)ut,ntraverse_p(ut)
!          write(*,*)ut,j,np,lospaths_p_up(ut)/lospaths(np)
!           write(*,14)'regular feb',(sunpaths(np,k),k=np,1,-1)
!           do j = 1, nfine
!             write(*,14)'regular feb',(sunpaths_fine(np,k,j),k=np,1,-1)
!             if (j.eq.partials_fineidx(ut))
!     &        write(*,14)'regular ut ',(sunpaths_p(ut,k),k=np,1,-1)
!           enddo
!           write(*,14)'regular feb',0.0d0,(sunpaths(np-1,k),k=np-1,1,-1)
!
!  Sun paths, Not direct sun , with tangent point. Code similar to above.
!  TANGR = tangent point radius for this ray
!   ntraverse_p(ut) = Number of layers traversed by ray.

!  Work downwards from TOA (sine rule) to level immediately above
!  the tangent layer. Don't forget to double the path length for
!  any layers which are traversed twice.

!  Tangent layer path length. Twice again. The following check is good.
!  check       write(*,*)tangr/dtan(th1),radii(krad)*dcos(th1)

          if (.not.direct_sun ) then
            tangr = stheta*radii(n)
            krad = nlayers
            do while (tangr.gt.radii(krad))
              krad = krad - 1
            enddo
            ntraverse_p(ut) = krad + 1
            sth0 = tangr/radii(0)
            th0 = dasin(sth0)
            do k = 1, krad
              sth1 = radii(k-1)*sth0/radii(k)
              th1 = dasin(sth1)
              ks1 = th1-th0
              fac = 1.0d0
              if ( k.gt.np) fac = 2.0d0
              sunpaths_p(ut,k) = fac*dsin(ks1)*radii(k-1)/sth1
              sth0 = sth1
              th0  = th1
            enddo
            sunpaths_p(ut,krad+1)=2.0d0*radii_p(ut)*dcos(th1)
          endif

!  Finish partials loop

        enddo
      endif

!  Finish

      return
END SUBROUTINE outgoing_sphergeom_fine_dn

!

SUBROUTINE multi_outgoing_adjustgeom                              &
         ( max_vza, max_sza, max_azm, n_vza, n_sza, n_azm,        & ! Input
          hsurface, eradius, do_adjust_surface,                   & ! Input
          alpha_boa, theta_boa, phi_boa,                          & ! Input
          alpha_ssa, theta_ssa, phi_ssa, fail, message, trace )     ! Output

! stand-alone geometry routine for adjusting the outgoing correction
!    starting inputs are the BOA values of SZA, VZA and PHI
!    need also the height of new surface, earth radius.

!  Height grid here is artificial

      implicit none

!  inputs

      integer       :: max_vza, max_sza, max_azm
      integer       :: n_vza, n_sza, n_azm
      real(fpk)     :: eradius, hsurface
      logical       :: do_adjust_surface
      real(fpk)     :: alpha_boa (max_vza)
      real(fpk)     :: theta_boa (max_sza)
      real(fpk)     :: phi_boa   (max_azm)

!  outputs

      real(fpk)     :: alpha_ssa (max_vza)
      real(fpk)     :: theta_ssa (max_vza,max_sza,max_azm)
      real(fpk)     :: phi_ssa   (max_vza,max_sza,max_azm)
      logical       :: fail
      character*(*) :: message, trace

!  Local

      integer       :: j, i, k
      real(fpk)     :: deg_to_rad, ex, ey, ez, px, py, pz
      real(fpk)     :: salpha_boa, calpha_boa, sphi_boa
      real(fpk)     :: stheta_boa, ctheta_boa, cphi_boa
      real(fpk)     :: ksi, cksi, sksi, xicum, cosscat_up
      real(fpk)     :: phi_all, alpha_all, theta_all
      real(fpk)     :: ctheta, stheta, calpha, salpha, cphi
      real(fpk)     :: b,rssa
      character*2   :: c2

!  Initialise output

      fail = .false.
      message = ' '
      trace   = ' '

!  check range of inputs

      do j = 1, n_vza
       if ( alpha_boa(j).ge.90.0d0.or.alpha_boa(j).lt.0.0d0 ) then
        write(c2,'(I2)')J
        message = 'boa LOS angle outside range [0,90])'
        trace   = 'Change Boa Los angle, number '//C2
        fail    = .true.
        return
       endif
      enddo

      do k = 1, n_azm
        if ( phi_boa(k).lt.0.0d0   ) phi_boa(k) = - phi_boa(k)
        if ( phi_boa(k).gt.180.0d0 ) phi_boa(k) = 360.0d0 - phi_boa(k)
      enddo

      do i = 1, n_sza
       if ( theta_boa(i).ge.90.0d0.or.theta_boa(i).lt.0.0d0 ) then
        write(c2,'(I2)')I
        message = 'Boa SZA angle outside range [0,90])'
        trace   = 'Change Boa SZA angle, number '//C2
        fail    = .true.
        return
       endif
      enddo

!  No adjustment, just copy and exit

      if ( .not. do_adjust_surface ) then
        do j = 1, n_vza
          alpha_ssa(j)   = alpha_boa(j)
          do i = 1, n_sza
            do k = 1, n_azm
              theta_ssa(j,i,k) = theta_boa(i)
              phi_ssa(j,i,k)   = phi_boa(k)
            enddo
          enddo
        enddo
        return
      endif

!  conversion

      deg_to_rad = dacos(-1.0d0) / 180.0d0

!  Radius of surface

      rssa = hsurface + eradius

!  Start VZA loop

      do j = 1, n_vza

        alpha_all = alpha_boa(j) * deg_to_rad
        salpha_boa = dsin(alpha_all)
        calpha_boa = dcos(alpha_all)

!  Special case. Direct nadir viewing. Compute everything and Exit.
!    (This is the same as the regular pseudo-spherical )

        if ( salpha_boa.eq.0.0d0 ) then
          alpha_ssa(j)   = alpha_boa(j)
          do i = 1, n_sza
            do k = 1, n_azm
              theta_ssa(j,i,k) = theta_boa(i)
              phi_ssa(j,i,k)   = phi_boa(k)
            enddo
          enddo
        endif

!  Los angle, general case

        if ( salpha_boa.ne.0.0d0 ) then

          salpha       = eradius * salpha_boa / rssa
          alpha_ssa(j) = dasin(salpha)
          calpha       = dcos(alpha_ssa(j))

!  Lospaths

          ksi  = alpha_all - alpha_ssa(j)
          sksi = dsin(ksi)
          cksi = dcos(ksi)
          xicum = ksi

!  output angle in degrees

          alpha_ssa(j) = alpha_ssa(j) / deg_to_rad

!  Los vector

          px = - rssa * dsin(xicum)
          py = 0.0d0
          pz =   rssa * dcos(xicum)

!  Start SZA loop

          do i = 1, n_sza

            theta_all = theta_boa(i) * deg_to_rad
            stheta_boa = dsin(theta_all)
            ctheta_boa = dcos(theta_all)

!  Start azimuth loop

            do k = 1, n_azm

              phi_all   = phi_boa(k)   * deg_to_rad
              cphi_boa  = dcos(phi_all)
              sphi_boa  = dsin(phi_all)

!  define Unit solar vector

              ex = - stheta_boa * cphi_boa
              ey = - stheta_boa * sphi_boa
              ez = - ctheta_boa

!  Sun angle

              b = ex*px + ey*py + ez*pz
              ctheta = -b/rssa
              stheta = dsqrt(1.0d0-ctheta*ctheta)
              theta_ssa(j,i,k) = dacos(ctheta)/deg_to_rad
              if ( ctheta.lt.0.0d0 ) then
                write(c2,'(I2)')J
                message = 'LOS-path SZA angle outside range [0,90])'
                trace   = 'Check inputs for LOS angle '//c2
                fail    = .true.
                return
              endif

!  scattering angle

              cosscat_up  = - calpha_boa * ctheta_boa + &
                              salpha_boa * stheta_boa * cphi_boa

!  Fix phi by using constancy of scatter angle

              if ( phi_boa(k).eq.180.0d0 ) then
                phi_ssa(j,i,k) = phi_all
              else if ( phi_boa(k) .eq. 0.0d0 ) then
                phi_ssa(j,i,k) = 0.0d0
              else
                cphi = (cosscat_up+calpha*ctheta)/stheta/salpha
                if ( cphi.gt.1.0d0) cphi = 1.0d0
                if ( cphi.lt.-1.0d0) cphi = -1.0d0
                phi_ssa(j,i,k) = dacos(cphi)
              endif
              phi_ssa(j,i,k) = phi_ssa(j,i,k) / deg_to_rad

!  End Azimuth and solar loops

            enddo
          enddo

!  End of regular LOS case

        endif

!  End los loop

      enddo

!  Finish

      return
END SUBROUTINE multi_outgoing_adjustgeom

!

SUBROUTINE ObsGeom_outgoing_adjustgeom                            &
         ( max_vza, max_sza, max_azm, n_vza, n_sza, n_azm,        & ! Input
          hsurface, eradius, do_adjust_surface,                   & ! Input
          alpha_boa, theta_boa, phi_boa,                          & ! Input
          alpha_ssa, theta_ssa, phi_ssa, fail, message, trace )     ! Output

! stand-alone geometry routine for adjusting the outgoing correction
!    starting inputs are the BOA values of SZA, VZA and PHI
!    need also the height of new surface, earth radius.

!  Height grid here is artificial

      implicit none

!  inputs

      integer       :: max_vza, max_sza, max_azm
      integer       :: n_vza, n_sza, n_azm
      real(fpk)     :: eradius, hsurface
      logical       :: do_adjust_surface
      real(fpk)     :: alpha_boa (max_vza)
      real(fpk)     :: theta_boa (max_sza)
      real(fpk)     :: phi_boa   (max_azm)

!  outputs

      real(fpk)     :: alpha_ssa (max_vza)
      real(fpk)     :: theta_ssa (max_vza,max_sza,max_azm)
      real(fpk)     :: phi_ssa   (max_vza,max_sza,max_azm)
      logical       :: fail
      character*(*) :: message, trace

!  Local

      integer       :: j, i, k
      real(fpk)     :: deg_to_rad, ex, ey, ez, px, py, pz
      real(fpk)     :: salpha_boa, calpha_boa, sphi_boa
      real(fpk)     :: stheta_boa, ctheta_boa, cphi_boa
      real(fpk)     :: ksi, cksi, sksi, xicum, cosscat_up
      real(fpk)     :: phi_all, alpha_all, theta_all
      real(fpk)     :: ctheta, stheta, calpha, salpha, cphi
      real(fpk)     :: b,rssa
      character*2   :: c2

!  Initialise output

      fail = .false.
      message = ' '
      trace   = ' '

!  check range of inputs

      do j = 1, n_vza
       if ( alpha_boa(j).ge.90.0d0.or.alpha_boa(j).lt.0.0d0 ) then
        write(c2,'(I2)')J
        message = 'boa LOS angle outside range [0,90])'
        trace   = 'Change Boa Los angle, number '//C2
        fail    = .true.
        return
       endif
      enddo

      do j = 1, n_azm
        if ( phi_boa(j).lt.0.0d0   ) phi_boa(j) = - phi_boa(j)
        if ( phi_boa(j).gt.180.0d0 ) phi_boa(j) = 360.0d0 - phi_boa(j)
      enddo

      do j = 1, n_sza
       if ( theta_boa(j).ge.90.0d0.or.theta_boa(j).lt.0.0d0 ) then
        write(c2,'(I2)')J
        message = 'Boa SZA angle outside range [0,90])'
        trace   = 'Change Boa SZA angle, number '//C2
        fail    = .true. ; return
       endif
      enddo

!  No adjustment, just copy and exit

      if ( .not. do_adjust_surface ) then
        do j = 1, n_vza
          alpha_ssa(j)     = alpha_boa(j)
          theta_ssa(1,j,1) = theta_boa(j)
          phi_ssa(1,1,j)   = phi_boa(j)
        enddo
        return
      endif

!  conversion

      deg_to_rad = dacos(-1.0d0) / 180.0d0

!  Radius of surface

      rssa = hsurface + eradius

!  Start Geometry loop (VZA index will be the driver.....)

      do j = 1, n_vza

        alpha_all = alpha_boa(j) * deg_to_rad
        salpha_boa = dsin(alpha_all)
        calpha_boa = dcos(alpha_all)

!  Special case. Direct nadir viewing. Compute everything and Exit.
!    (This is the same as the regular pseudo-spherical )

        if ( salpha_boa.eq.0.0d0 ) then
          alpha_ssa(j)     = alpha_boa(j)
          theta_ssa(1,j,1) = theta_boa(j)
          phi_ssa(1,1,j)   = phi_boa(j)
        endif

!  Los angle, general case

        if ( salpha_boa.ne.0.0d0 ) then

          salpha       = eradius * salpha_boa / rssa
          alpha_ssa(j) = dasin(salpha)
          calpha       = dcos(alpha_ssa(j))

!  Lospaths

          ksi  = alpha_all - alpha_ssa(j)
          sksi = dsin(ksi)
          cksi = dcos(ksi)
          xicum = ksi

!  output angle in degrees

          alpha_ssa(j) = alpha_ssa(j) / deg_to_rad

!  Los vector

          px = - rssa * dsin(xicum)
          py = 0.0d0
          pz =   rssa * dcos(xicum)

!  SZA calculation

         theta_all = theta_boa(j) * deg_to_rad
         stheta_boa = dsin(theta_all)
         ctheta_boa = dcos(theta_all)

!  azimuth calculation

         phi_all   = phi_boa(j)   * deg_to_rad
         cphi_boa  = dcos(phi_all)
         sphi_boa  = dsin(phi_all)

!  define Unit solar vector

         ex = - stheta_boa * cphi_boa
         ey = - stheta_boa * sphi_boa
         ez = - ctheta_boa

!  Sun angle

         b = ex*px + ey*py + ez*pz
         ctheta = -b/rssa
         stheta = dsqrt(1.0d0-ctheta*ctheta)
         theta_ssa(1,j,1) = dacos(ctheta)/deg_to_rad
         if ( ctheta.lt.0.0d0 ) then
            write(c2,'(I2)')J
            message = 'LOS-path SZA angle outside range [0,90])'
            trace   = 'Check inputs for LOS angle '//c2
            fail    = .true. ; return
         endif

!  scattering angle

         cosscat_up  = - calpha_boa * ctheta_boa + &
                         salpha_boa * stheta_boa * cphi_boa

!  Fix phi by using constancy of scatter angle

         if ( phi_boa(j).eq.180.0d0 ) then
            phi_ssa(1,1,j) = phi_all
         else if ( phi_boa(j) .eq. 0.0d0 ) then
            phi_ssa(1,1,j) = 0.0d0
         else
            cphi = (cosscat_up+calpha*ctheta)/stheta/salpha
            if ( cphi.gt.1.0d0) cphi = 1.0d0
            if ( cphi.lt.-1.0d0) cphi = -1.0d0
            phi_ssa(1,1,j) = dacos(cphi)
         endif
         phi_ssa(1,1,j) = phi_ssa(1,1,j) / deg_to_rad

!  End of regular LOS case

        endif

!  End los loop

      enddo

!  Finish

      return
END SUBROUTINE ObsGeom_outgoing_adjustgeom

!

SUBROUTINE losonly_outgoing_adjustgeom &
       ( max_vza, n_vza, hsurface, eradius, do_adjust_surface, alpha_boa, &
         alpha_ssa, fail, message, trace )

      implicit none

!  Thermal only (no solar)

! Stand-alone geometry routine for adjusting the outgoing correction
!    starting inputs are the boa values of vza
!    need also the height of new surface, earth radius.

!  Height grid here is artificial

!  Inputs

      integer, intent(in) ::           max_vza
      integer, intent(in) ::           n_vza
      real(fpk), intent(in) ::         eradius, hsurface
      logical, intent(in) ::           do_adjust_surface
      real(fpk), intent(in) ::         alpha_boa (max_vza)

!  Outputs

      real(fpk), intent(out) ::        alpha_ssa (max_vza)
      logical, intent(out) ::          fail
      character (len=*), intent(inout) :: message, trace

!  Local

      integer ::    j
      real(fpk) ::  deg_to_rad, pie, pi2
      real(fpk) ::  salpha_boa, calpha_boa
      real(fpk) ::  alpha_all,  calpha, salpha, rssa
      character (len=2) :: c2

!  Initialise output

      fail    = .false.
      message = ' '
      trace   = ' '

!  Check range of inputs

      do j = 1, n_vza
       if ( alpha_boa(j).ge.90.0d0.or.alpha_boa(j).lt.0.0d0 ) then
        write(c2,'(i2)')j
        message = 'boa los angle outside range [0,90])'
        trace   = 'change boa los angle, number '//c2
        fail    = .true.
        return
       endif
      enddo

!  No adjustment, just copy and exit

      if ( .not. do_adjust_surface ) then
        do j = 1, n_vza
          alpha_ssa(j)   = alpha_boa(j)
        enddo
        return
      endif

!  Conversion

      pie = dacos(-1.0d0)
      pi2 = 2.0d0 * pie
      deg_to_rad = pie / 180.0d0

!  Radius of surface

      rssa = hsurface + eradius

!  Start vza loop

      do j = 1, n_vza

        alpha_all = alpha_boa(j) * deg_to_rad
        salpha_boa = dsin(alpha_all)
        calpha_boa = dcos(alpha_all)

!  Special case. direct nadir viewing. compute everything and exit.
!    (this is the same as the regular pseudo-spherical )

        if ( salpha_boa.eq.0.0d0 ) then
          alpha_ssa(j)   = alpha_boa(j)
        else
           salpha       = eradius * salpha_boa / rssa
           alpha_ssa(j) = dasin(salpha)
           calpha       = dcos(alpha_ssa(j))
           alpha_ssa(j) = alpha_ssa(j) / deg_to_rad
        endif

!  End los loop

      enddo

!  Finish

      return
END SUBROUTINE losonly_outgoing_adjustgeom

!  Finish

end module lidort_geometry

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


MODULE GRTOPH

!! Module for transformation from phase space into physical space

  USE PHTOGR

CONTAINS

  SUBROUTINE GRPH213(CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)

! WRONG>>> DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
! RAUM AUF KUGELKOORDINATEN DURCH

! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
! CXM   = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
! WSAVE  = Working Array fuer Fouriertransformation
! Z   = LEGENDREFUNKTIONSWERTE
!
! MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
! MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
! MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
! MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN

    IMPLICIT REAL (A-H,O-Z)

!   Anzahl der Gitterpunkte pro Breitenkreis des reduzierten
!   Gauss'schen Gitters
    INTEGER MLAT(MAXB),ISIZE,IFAX(10,MAXB)

!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
    REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)

!   LOGICAL*1 USED(((216*217)/2+1)*160)

    DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
    REAL FELD(MAXL,MLEVEL)
    DIMENSION WSAVE(8*MAXB+15,MAXB/2)
    REAL W(MAXB)
    DIMENSION IND(MAXB)

    IND(1)=0
    DO 6 J=2,MAXB/2
      IND(j)=IND(J-1)+MLAT(J-1)
6   CONTINUE
!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
    DO 16 L=1,MLEVEL
      CALL GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
16  CONTINUE
!$omp end parallel do

    RETURN

  END SUBROUTINE GRPH213

  SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)

!! Carries out the transformation from the phase space into the physical
!! space on spherical coordinates

! CXMN  = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
!         CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
! CXM   = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
! FELD  = FELD DER METEOROLOGISCHEN VARIABLEN
! WSAVE = Working Array fuer Fouriertransformation
! Z     = LEGENDREFUNKTIONSWERTE
!
! MNAUF ANZAHL DER FOURIERKOEFFIZIENTEN
! MAXL  ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
! MAXB  ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN

    IMPLICIT REAL (A-H,O-Z)

!   FELD DER FOURIERKOEFFIZIENTEN
    REAL CXMS(4*(MNAUF+1))
    REAL CXMA(4*(MNAUF+1))
    REAL,ALLOCATABLE :: CXM(:,:)

!   Anzahl der Gitterpunkte pro Breitenkreis des reduzierten
!   Gauss'schen Gitters
    INTEGER MLAT(MAXB),ISIZE

!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
    REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)

!   LOGICAL*1 USED(((216*217)/2+1)*160)

    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
    REAL FELD(MAXL,MLEVEL)
    REAL WSAVE(8*MAXB+15,MAXB/2)
    INTEGER IFAX(10,MAXB)
    REAL W(MAXB)
    INTEGER IND(MAXB)

    ALLOCATE(CXM( 4*MAXB,MAXB))
    DO 5 J=1,MAXB/2
      CXMS(1:MLAT(J))=FELD(IND(J)+1:IND(J)+MLAT(J),L)
      CALL RFOUFTR(CXMS,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
      CXMA(1:MLAT(J))=FELD(MAXL-IND(J)-MLAT(J)+1:MAXL-IND(J),L)
      CALL RFOUFTR(CXMA,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
      DO 4 I=1,2*(MNAUF+1)
        CXM(I,J)=CXMS(I)+CXMA(I)
        CXM(I,MAXB+1-J)=CXMS(I)-CXMA(I)
4     CONTINUE
5   CONTINUE
    CALL LGTR213(CXMN(0,L),CXM,Z,W,MLAT,MNAUF,MAXB)

    DEALLOCATE(CXM)

    RETURN

  END SUBROUTINE GRPHSUB
!
  SUBROUTINE LGTR213(CXMN,CXM,Z,W,MLAT,MNAUF,MAXB)

!! Calculates the KFFKs (Fourier coefficients) CXMN

    IMPLICIT REAL (A-H,O-Z)
    INTEGER MLAT(MAXB)
    DIMENSION CXM(0:4*MAXB-1,MAXB)
    DIMENSION CXMN(0:2*(((MNAUF+1)*MNAUF)/2+MNAUF)+1)
    REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
    REAL*8 W(MAXB),CR,CI,HILF
    LOGICAL EVEN

    LL=0
    LLP=0
    DO 1 I=0,MNAUF
      KM=0
9     KM=KM+1
      IF (MLAT(KM) .LE. 2*I) THEN
        GOTO 9
      END IF
      DO 2 J=I,MNAUF
        CR=0
        CI=0
        EVEN=MOD(I+J,2) .EQ. 0
        IF (EVEN) THEN
          DO 3 K=KM,MAXB/2
            HILF=W(K)*Z(K,LLP)
            CR=CR+CXM(2*I,K)*HILF
            CI=CI+CXM(2*I+1,K)*HILF
3         CONTINUE
        ELSE
          DO 4 K=KM,MAXB/2
            HILF=W(K)*Z(K,LLP)
            CR=CR+CXM(2*I,MAXB+1-K)*HILF
            CI=CI+CXM(2*I+1,MAXB+1-K)*HILF
4         CONTINUE
        END IF
5       CXMN(2*LL)=CR
        CXMN(2*LL+1)=CI
        LL=LL+1
        LLP=LLP+1
2     CONTINUE
      LLP=LLP+2
1   CONTINUE
    RETURN

  END SUBROUTINE LGTR213

  SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)

!! Calculates the Fourier sum with an FFT algorithm

    IMPLICIT REAL (A-H,O-Z)
    DIMENSION CXM(0:2*MAXL-1)
    DIMENSION FELD(MAXL),TRIGS(2*MAXL)
    DIMENSION WSAVE(MAXAUF)
    INTEGER IFAX(10)

! NORMIERUNG...
    WSAVE(1)=CXM(MAXL-1)

    CXM(1:MAXL)=CXM(0:MAXL-1)/2
    CXM(0)=WSAVE(1)/2
!    CALL CFFTF(MAXL,CXM,WSAVE)
    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,-1)
    RETURN
  END SUBROUTINE RFOUFTR

END MODULE GRTOPH