**************************************************************************
***   pcaknmi.f == 4 April 2005
***
***   polarisation correction subroutines of gomecal -- version 1.2
***
***   The routines originate from Nick Schutgens, KNMI,
***   and are adapted for use with gomecal.
***
***   contents: subroutines CorrectAngles, Validity, PCA_CorrFactor
***                         akima_ex, akima_su
***             functions   Parametrise_SS, Parametrise_SG, GDF, akima_po
**************************************************************************
*
*  The original routines were written in Fortran-90. These have been
*  converted to Fortran-77 to harmonise them with the rest of the
*  gomecal code. The definitions in the original module "Constants"
*  are put in the routines where they are used.
*
*  The following lines are the header from the original routines; the
*  original individual routine descriptions have been copied too.
*
!==========================================================================
!      Module PCA_KNMI v1.2
!
!  This module contains all functions and subroutines necessary for
!  the improved polarisation correction of GOME level 0 data.
!  As it uses parametrisations for UV polarisation, it performs
!  substantially better than the operational GOME lvl 1 data processor
!  (up to v2.0 at least) in the UV (290-330 nm).
!
!  The user only needs to call <PCA_CorrFactor>
!  and supply it with:
!   - cosine solar zenith angle for the center of a GOME pixel (B point)
!   - cosine viewing zenith angle for the center of a GOME pixel (B point)
!   - albedo, either surface albedo or (in case of clouds) a scene albedo
!   - total ozone column (in Dobson units)
!   - polarisation measurements from GOME and the theoretical single
!        Rayleigh scattering polarisation
!
!  The paper by N.A.J. Schutgens & P. Stammes referred to is:
!  'Parametrisation of Earth's polarisation spectrum from 290 to 330 nm'
!  J. Quant. Spectrosc. & Radiat. Transfer, vol. 75, pp. 239--255, 2002.
!
!  Created: 15-MAY-2001
!  Author: Nick Schutgens, KNMI
!--------------------------------------------------------------------------
!  24-JULY-2001: - Corrected limits on DO loops in akima_ex/akima_su
!                  causing array boundary overflows
!  23-JULY-2001: - Changed REAL to REAL(KIND=8)
!                - Changed FUNCTION p_CorrFactor to SUBROUTINE PCA_CorrFactor
!                - Added SUBROUTINE CorrectAngles
!==========================================================================
*
**************************************************************************
*
      subroutine CorrectAngles(SolSZA,ViewSZA,RelAzi,SatH,EarthR)
c
!==========================================================================
!  Subroutine :  CorrectAngles
!  Input      :  SolSZA   = solar zenith angle B w.r.t. North -- [1]
!                ViewSZA  = view zenith angle B w.r.t. satellite -- [1]
!                RelAzi   = relative azimuth B ( =< 180 degrees) -- [1]
!                SatH     = satellite height [km] -- [1]
!                EarthR   = Earth curvature radius [km] -- [1]
!  Output     :  SolSZA   = solar zenith angle, at local tangent plane
!                ViewSZA  = view zenith angle, at local tangent plane
!                RelAzi   = relative azimuth angle, at local tangent plane
!         [1] = see GOME product
!  Description:  corrects scattering geometry of GOME product to angles for
!                local tangent plane
!  References :  none
!  Libraries  :  none
!  Created    :  23-JULY-2001
!  Version nr :  1.0
!  Author     :  Nick Schutgens, KNMI
!  Comments   :  adaptations by Jos van Geffen:
!                 - write "ViewSZA=ViewSZA*pi/180.0" (etc.) instead of the
!                   not uniquely defined ViewSZA=pi/180.0*ViewSZA"
!                 - remove the "pause" and instead write a message to the
!                   logfile (cf. subroutine Validate).
!--------------------------------------------------------------------------
c
      implicit none
c
c  Common blocks.
c
      include 'gomecal.inc'
c
c  Local variables.
c
      real*8 SolSZA,ViewSZA,RelAzi,SatH,EarthR
      real*8 rpi,RotAngle,cosSolSZA
      parameter(rpi=3.141592653589793d0)
c
c ------------------
c
c  Redefine all angles in radians
c
      ViewSZA=ViewSZA*rpi/180.d0
      SolSZA =SolSZA*rpi/180.d0
      RelAzi =RelAzi*rpi/180.d0
c
c  First obtain rotation angle
c
      RotAngle=-ViewSZA+dASIN( (EarthR+SatH)/EarthR*dSIN(ViewSZA) )
c
      if ((RotAngle.lt.0.d0).and.(llog)) write(77,10)RotAngle
  10  format('      --> rotation angle RotAngle is negative:',f5.2,
     #'; continuing')
c
c  Correct the solar SZA angle
c
      cosSolSZA=dCOS(SolSZA)*dCOS(RotAngle)+
     #          dSIN(SolSZA)*dSIN(RotAngle)*dCOS(RelAzi)
c
c  Correct relative azimuth angle
c
      if (dSQRT(1.d0-cosSolSZA**2).gt.1.d-4) then
        RelAzi=dASIN(dSIN(SolSZA)*dSIN(RelAzi)/dSQRT(1.d0-cosSolSZA**2))
       else
        RelAzi=RelAzi
      endif
c
c  Correct viewing angle
c
      ViewSZA=dASIN( (EarthR+SatH)/EarthR*dSIN(ViewSZA) )
c
c  Redefine all angles in degrees
c
      ViewSZA=ViewSZA*180.d0/rpi
      SolSZA =dACOS(cosSolSZA)*180.d0/rpi
      RelAzi =RelAzi*180.d0/rpi
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      subroutine Validity(Mu0,Mu,Albedo,O3Column)
c
!==========================================================================
!  SUBROUTINE : Validity
!  Input      : Mu0       (cosine solar zenith angle)
!               Mu        (cosine viewing zenith angle)
!               Albedo    (surface albedo)
!               O3Column  (total ozone column)
!  Output     : none
!  Description: checks if the input variables are within the bounds for which
!               the parametrisation was developed. Out of bound values do not
!               necessarily imply a bad parametrisation. Use common sense.
!  Referrences: see paper by NAJ Schutgens & P Stammes
!  Libraries  : none
!  Created    : 19-JULY-2001
!  Version nr : 1.0
!  Author     : NAJ Schutgens
!  Comments   : none
!             : layout of output adapted for use with gomecal, which calls
!               this routine only if the logical llog is .true., meaning
!               that the log file is connected to unit 77. --  JvG
!--------------------------------------------------------------------------
c
      implicit none
c
c  Common blocks.
c
      include 'gomecal.inc'
c
c  Local variables
c
      real*8 Mu0,Mu,Albedo,O3Column
      real*8 Mu0_min,Mu0_max,Mu_min,Mu_max
      real*8 Albedo_min,Albedo_max,O3Column_min,O3Column_max
c
c ------------------
c
      Mu0_min=0.25d0
      Mu0_max=1.0d0
      if ((Mu0.lt.Mu0_min).or.(Mu0.gt.Mu0_max)) then
         if (llog) write(77,11)Mu0_min,Mu0_max,Mu0
      endif
c
      Mu_min=0.8d0
      Mu_max=1.0d0
      if ((Mu.lt.Mu_min).or.(Mu.gt.Mu_max)) then
         if (llog) write(77,12)Mu_min,Mu_max,Mu
      endif
c
      Albedo_min=0.0d0
      Albedo_max=1.0d0
      if ((Albedo.lt.Albedo_min).or.(Albedo.gt.Albedo_max)) then
         if (llog) write(77,13)Albedo_min,Albedo_max  !,Albedo
      endif
c
      O3Column_min=222.0d0
      O3Column_max=439.0d0
      if ((O3Column.lt.O3Column_min).or.(O3Column.gt.O3Column_max)) then
         if (llog) write(77,14)O3Column_min,O3Column_max  !,O3Column
      endif
c
  11  format('      --> parametrisation derived for ',f5.3,
     #' <= mu0 <= ',f5.3,'; file has ',f5.3)
  12  format('      --> parametrisation derived for ',f5.3,
     #' <= mu  <= ',f5.3,'; file has ',f5.3)
  13  format('      --> parametrisation derived for ',f5.3,
     #' <= albedo   <= ',f5.3)                   !,'; file has ',f5.3)
  14  format('      --> parametrisation derived for ',f5.1,
     #' <= O3column <= ',f5.1)                   !,'; file has ',f5.1)
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      subroutine PCA_CorrFactor(Mu0,Mu,Albedo,O3Column,gome_l,gome_p,
     #                          Lambda,Eta,nLambda,CorrFac)
c
!==========================================================================
!  Function   : PCA_CorrFactor
!  Input      : Mu0      = cosine solar zenith angle
!               Mu       = cosine viewing zenith angle
!               Albedo   = surface albedo
!               O3Column = total ozone column
!               gome_l   = wavelengths of GOME polarisation data
!               gome_p   = fractional polarisation values of GOME data
!               Lambda   = wavelength grid for which interpolation is sought
!               Eta      = polarisation sensitivity, corrected for scan
!                          mirror angle
!               nLambda  = nr of wavelengths
!  Output     : CorrFac  = correction factor for polarisation sensitivity
!  Description: using GOME polarisation data, an interpolated fractional
!               polarisation is found from which a correction factor for
!               the polarisation sensitivity is found
!  References : see paper by NAJ Schutgens & P Stammes and Balzer et al.
!               (see GOME level 0 to 1 Algorithms description,
!               ER-TN_DLR_GO_0022, DLR, 1996)
!  Libraries  : none
!  Created    : 15-MAY-2001
!  Version nr : 1.0
!  Author     : NAJ Schutgens
!  Comments   : This algorithm is different from the one present in
!               the GOME extractor software: 1) it starts with the GDF and
!               then matches the Akima interpolation unto it 2) it uses a
!               different strategy in case of missing polarisation data
!             : variable Lambda0 was output-ed first too, but that value is
!               not needed anyfurther, hence it is removed -- JvG
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      integer nLambda
      real*8 Mu0,Mu,Albedo,O3Column
      real*8 gome_l(8),gome_p(8)
      real*8 Lambda(nLambda),Eta(nLambda),CorrFac(nLambda)
c
      integer theor_VIS,theor_UV,max_n_Akima,PMD1,PMD3
      parameter(theor_VIS=7,theor_UV=8,max_n_Akima=10,PMD1=1,PMD3=3)
c -- PMD2=2 was defined in original, but is not used, so omitted here
      real*8 Lambda_EndGDF,Lambda_Subst,dLambda_s,dLambda_l
      parameter(Lambda_EndGDF=325.d0,Lambda_Subst=700.d0)
      parameter(dLambda_s=5.d0,dLambda_l=20.d0)
      real*8 Parametrise_SS,Parametrise_SG,GDF,akima_po
c
      real*8 x_Akima(max_n_Akima),y_Akima(max_n_Akima)
      real*8 a_Akima(max_n_Akima),b_Akima(max_n_Akima)
      real*8 c_Akima(max_n_Akima),d_Akima(max_n_Akima)
      real*8 w0,DeltaLambda,p_00,Lambda0
      integer n_Akima,i
c
c ------------------
c
c  Perform parametrisation.
c
      Lambda0 = Parametrise_SS(Mu0,Mu,Albedo,O3Column)
      DeltaLambda = - ( Parametrise_SG(Mu0,Mu,Albedo,O3Column)-Lambda0 )
     #                / dLOG(2.d0-dSQRT(3.d0))
c
c  Determine valid fractional polarisation values and create substitute
c  if necessary.
c
      n_Akima=0
c
c  Only if a valid theoretical p (VIS integration time) is present
c  do we check the PMD values.
c
      if ((gome_p(theor_VIS).ge.0.d0) .and.
     #    (gome_p(theor_VIS).le.1.d0)) then
c
         do i=PMD1,PMD3
            if ((gome_p(i).ge.0.d0) .and. (gome_p(i).le.1.d0)) then
               n_Akima=n_Akima+1
               x_Akima(n_Akima)=gome_l(i)
               y_Akima(n_Akima)=gome_p(i)
            endif
         enddo
c
c  If PMD values present, use first PMD to determine offset of GDF.
c
         if (n_Akima.gt.0) then
            p_00=y_Akima(1)
            w0=4.d0*(gome_p(theor_VIS)-p_00)
c
c  If no PMD values present, make an educated guess for the offset
c  of the GDF.
c
          else
            p_00=(gome_p(theor_VIS)+0.5d0)/2.d0
            w0=4.d0*(gome_p(theor_VIS)-p_00)
            n_Akima=1
            x_Akima(n_Akima)=Lambda_Subst
            y_Akima(n_Akima)=p_00
         endif
c
c  If only theoretical p (UV integration time) present, make an educated
c  guess for the offset of the GDF.
c
       else if ((gome_p(theor_UV).ge.0.d0) .and.
     #          (gome_p(theor_UV).le.1.d0)) then
c
         p_00=(gome_p(theor_UV)+0.5d0)/2.d0
         w0=4.d0*(gome_p(theor_UV)-p_00)
         n_Akima=1
         x_Akima(n_Akima)=Lambda_Subst
         y_Akima(n_Akima)=p_00
c
c  If no theoretical p values present, resort to damage control.
c
       else
         p_00=0.5d0
         w0=0.d0
         n_Akima=1
         x_Akima(n_Akima)=Lambda_Subst
         y_Akima(n_Akima)=p_00
c
      endif
c
c  Extend Akima array left and right with horizontal branches
c
      do i=n_Akima,1,-1
         x_Akima(i+3)=x_Akima(i)
         y_Akima(i+3)=y_Akima(i)
      enddo
      do i=1,3
         x_Akima(i)=Lambda_EndGDF-(2-i)*dLambda_s
         y_Akima(i)=GDF(x_Akima(i),Lambda0,DeltaLambda,w0,p_00)
      enddo
      do i=n_Akima+4,n_Akima+5
         x_Akima(i)=x_Akima(n_Akima+3)+(i-n_Akima+3)*dLambda_l
         y_Akima(i)=y_Akima(n_Akima+3)
      enddo
      n_Akima=n_Akima+5
c
c  Determine Akima interpolation coefficients.
c
      CALL akima_su(x_Akima,y_Akima,n_Akima,a_Akima,
     #              b_Akima,c_Akima,d_Akima)
c
c  Determine interpolated fractional polarisation.
c
      do i=1,nLambda
         if (Lambda(i).le.Lambda_EndGDF) then
            CorrFac(i)=GDF(Lambda(i),Lambda0,DeltaLambda,w0,p_00)
          else
            CorrFac(i)=akima_po(n_Akima,x_Akima,a_Akima,b_Akima,
     #                          c_Akima,d_Akima,Lambda(i))
         endif
      enddo
c
c  Determine polarisation correction factor.
c
      do i=1,nLambda
         CorrFac(i)=0.5d0*(1.d0+Eta(i)) /
     #              (CorrFac(i)*(1.d0-Eta(i))+Eta(i))
      enddo
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      real*8 function Parametrise_SS(Mu0,Mu,Albedo,O3Column)
c
!==========================================================================
!  Function   :  Parametrise_SS
!  Input      :  Mu0       (cosine solar zenith angle)
!                Mu        (cosine viewing zenith angle)
!                Albedo    (surface albedo)
!                O3Column  (total ozone column [DU])
!  Output     :  a single wavelength [nm]
!  Description:  the largest wavelength for which single Rayleigh scattering
!                still yields a good approximation of the fractional
!                polarisation, is parametrised in the four input variables
!  References :  see paper by NAJ Schutgens & P Stammes
!  Libraries  :  none
!  Created    :  15-MAY-2001
!  Version nr :  1.0
!  Author     :  Nick Schutgens, KNMI
!  Comments   :  In case of cloudy scenes, surface albedo may be replaced
!                by scene albedo for better results
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      real*8 Mu0,Mu,Albedo,O3Column
      real*8 Geom_Airmass,O3_fraction,O3ColumnRef
      real*8 C_SS_MA(0:2,0:2),C_SS_O3(0:2)
      integer i,j
c
c  Original parameter values.
c
c      data C_SS_MA / 308.6d0,-28.7d0,10.9d0,  !  C_SS_MA(i,0)
c     #                 0.0d0,  0.0d0, 0.0d0,  !  C_SS_MA(i,1)
c     #                 0.0d0,  0.0d0, 0.0d0 / !  C_SS_MA(i,2)
c      data C_SS_O3 / 0.0d0,7.6d0,-4.3d0 /
c
c  Updated parameter values, in accordance with the published paper.
c
      data C_SS_MA / 308.68d0,-29.10d0,11.46d0,  !  C_SS_MA(i,0)
     #                 0.00d0,  0.00d0, 0.00d0,  !  C_SS_MA(i,1)
     #                 0.00d0,  0.00d0, 0.00d0 / !  C_SS_MA(i,2)
      data C_SS_O3 / 0.00d0,7.58d0,-4.26d0 /
c
c ------------------
c
      O3ColumnRef=345.8d0
      Geom_Airmass=(Mu0+Mu)/(Mu0*Mu)
      O3_fraction=(O3Column-O3ColumnRef)/O3ColumnRef
      Parametrise_SS=0.d0
c
      do j=0,2
         do i=0,2
            Parametrise_SS = Parametrise_SS +
     #                     C_SS_MA(i,j)*(1.d0/Geom_Airmass)**i*Albedo**j
         enddo
         Parametrise_SS = Parametrise_SS + C_SS_O3(j)*O3_fraction**j
      enddo
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      real*8 function Parametrise_SG(Mu0,Mu,Albedo,O3Column)
c
!==========================================================================
!  Function   :  Parametrise_SG
!  Input      :  Mu0       (cosine solar zenith angle)
!                Mu        (cosine viewing zenith angle)
!                Albedo    (surface albedo)
!                O3Column  (total ozone column [DU])
!  Output     :  a single wavelength [nm]
!  Description:  the wavelength for which the global shape of the fractional
!                polarisation in the UV (as defined by the Generalised
!                Distribution Function) has its largest derivative,
!                is parametrised in the four input variables
!  References :  see paper by NAJ Schutgens & P Stammes
!  Libraries  :  none
!  Created    :  15-MAY-2001
!  Version nr :  1.0
!  Author     :  Nick Schutgens, KNMI
!  Comments   :  In case of cloudy scenes, surface albedo may be replaced
!                by scene albedo for better results
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      real*8 Mu0,Mu,Albedo,O3Column
      real*8 Geom_Airmass,O3_fraction,O3ColumnRef
      real*8 C_SG_MA(0:2,0:2),C_SG_O3(0:2)
      integer i,j
c
c  Original parameter values.
c
c      data C_SG_MA / 317.1d0,-46.2d0, 35.8d0,  !  C_SG_MA(i,0)
c     #                -2.1d0, 13.8d0,-18.2d0,  !  C_SG_MA(i,1)
c     #                 0.7d0, -9.3d0,  9.3d0 / !  C_SG_MA(i,2)
c      data C_SG_O3 / 0.0d0,7.277d0,-4.085d0 /
c
c  Updated parameter values, in accordance with the published paper.
c
      data C_SG_MA / 316.43d0,-41.89d0,29.49d0,  !  C_SG_MA(i,0)
     #                 0.33d0, -0.06d0, 0.66d0,  !  C_SG_MA(i,1)
     #                -1.11d0,  0.56d0,-3.46d0 / !  C_SG_MA(i,2)
      data C_SG_O3 / 0.00d0,7.20d0,-4.08d0 /
c
c ------------------
c
      O3ColumnRef=345.8d0
      Geom_Airmass=(Mu0+Mu)/(Mu0*Mu)
      O3_fraction=(O3Column-O3ColumnRef)/O3ColumnRef
      Parametrise_SG=0.d0
c
      do j=0,2
         do i=0,2
            Parametrise_SG = Parametrise_SG+
     #                     C_SG_MA(i,j)*(1.d0/Geom_Airmass)**i*Albedo**j
         enddo
         Parametrise_SG = Parametrise_SG + C_SG_O3(j)*O3_fraction**j
      enddo
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      real*8 function GDF(Lambda,Lambda0,DeltaLambda,w0,p_00)
c
!==========================================================================
!  Function   : GDF (Generalised distribution Function)
!  Input      : Lambda        wavelength
!               Lambda0       largest single scattering wavelength
!               DeltaLambda   defines steepness of GDF
!               w0            amplitude of GDF
!               p_00          offset of GDF, w0/4+p_00 is UV (290 nm)
!                             polarisation
!  Output     : value of fractional polarisation
!  Description: the GDF accurately describes the fractional polarisation
!               in the UV (290-330 nm)
!  References : see paper by NAJ Schutgens & P Stammes
!  Libraries  : none
!  Created    : 15-MAY-2001
!  Version nr : 1.0
!  Author     : Nick Schutgens, KNMI
!  Comments   : after an idea by Balzer et al. (see GOME level 0 to 1
!               Algorithms description, ER-TN_DLR_GO_0022, DLR, 1996)
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      real*8 Lambda,Lambda0,DeltaLambda,w0,p_00
c
c ------------------
c
      if (Lambda.le.Lambda0) then
         GDF=w0/4.d0+p_00
       else
         GDF=w0*dEXP(-(Lambda-Lambda0)/DeltaLambda) /
     #       (1.d0+dEXP(-(Lambda-Lambda0)/DeltaLambda))**2 + p_00
      endif
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      subroutine akima_ex(n_source,x_arr,y_arr,st)
c
!==========================================================================
!  Subroutine : akima_ex
!  Input      : n_source  nr of function values
!               x_arr     abscissa values
!               y_arr     function values
!  Output     : slopes of the curve
!  Description: extends a tabulated function that will be interpolated two
!               points left and right and calculates derivatives
!  References : GOME level 1 extractor source code
!  Libraries  : none
!  Created    : 15-MAY-2001
!  Version nr : 1.0
!  Author     : Nick Schutgens, KNMI
!  Comments   : Originally by Akima. C implementation by Bernd Aberle,
!               DLR/WT-DA-SE, from which this F90 code was developed
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      integer n_source
      real*8 x_arr(n_source+4),y_arr(n_source+4),st(n_source+4)
      real*8 tmp
      integer i
c
c ------------------
c
c  Auxiliary points are moved 2 indices up.
c
      do i=0,n_source-1
         x_arr((n_source+2)-i) = x_arr(n_source-i)
         y_arr((n_source+2)-i) = y_arr(n_source-i)
      enddo
c
c  Calculation of extrapolated x-values left and right.
c
      x_arr(2) = x_arr(3) + x_arr(4) - x_arr(5)
      x_arr(1) = x_arr(2) + x_arr(3) - x_arr(4)
      x_arr(n_source+3) = x_arr(n_source+2) + x_arr(n_source+1)
     #                    - x_arr(n_source)
      x_arr(n_source+4) = x_arr(n_source+3) + x_arr(n_source+2)
     #                    - x_arr(n_source+1)
c
c  Calculation of slopes of curves.
c
      do i=3,n_source+1
         tmp = x_arr(i+1) - x_arr(i)
         if (tmp.lt.1.d-10) then
            st(i) = 0.d0
          else
            st(i) = (y_arr(i+1) - y_arr(i)) / tmp
         endif
      enddo
c
c  Calculation of extrapolated y-values left and right.
c
      tmp = x_arr(3) - x_arr(2)
      y_arr(2) = tmp * (st(4) - 2.d0*st(3)) + y_arr(3)
      if ( tmp.lt.1.d-10) then
         st(2) = 0.d0
       else
         st(2) = (y_arr(3) - y_arr(2)) / tmp
      endif
c
      tmp = x_arr(2) - x_arr(1)
      y_arr(1) = tmp * (st(3) - 2.d0*st(2)) + y_arr(2)
      if (tmp.lt.1.d-10) then
         st(1) = 0.d0
       else
         st(1) = (y_arr(2) - y_arr(1)) / tmp
      endif
c
      tmp = x_arr(n_source+3) - x_arr(n_source+2)
      y_arr(n_source+3) = tmp * (2.d0*st(n_source+1) - st(n_source))
     #                    + y_arr(n_source+2)
      if (tmp.lt.1.d-10) then
         st(n_source+2) = 0.d0
       else
         st(n_source+2) = (y_arr(n_source+3) - y_arr(n_source+2)) / tmp
      endif
c
      tmp = x_arr(n_source+4) - x_arr(n_source+3)
      y_arr(n_source+4) = tmp * (2.d0*st(n_source+2) - st(n_source+1))
     #                    + y_arr(n_source+3)
      if (tmp.lt.1.d-10) then
         st(n_source+3) = 0.d0
       else
         st(n_source+3) = (y_arr(n_source+4) - y_arr(n_source+3)) / tmp
      endif
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      subroutine akima_su(x_arr,y_arr,n_source,a_arr,b_arr,c_arr,d_arr)
c
!==========================================================================
!  Subroutine : akima_su
!  Input      : n_source  nr of function values
!               x_arr     abscissa values
!               y_arr     function values
!  Output     : a_arr, b_arr, c_arr, d_arr  polynomial constants that define
!               Akima interpolation
!  Description: Determines the parameters of the Akima interpolation for a
!               tabulated function
!  References : GOME level 1 extractor source code
!  Libraries  : none
!  Created    : 15-MAY-2001
!  Version nr : 1.0
!  Author     : Nick Schutgens, KNMI
!  Comments   : Originally by Akima. C implementation by Bernd Aberle,
!               DLR/WT-DA-SE from which this F90 code was developed
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      integer n_source
      real*8 x_arr(n_source),y_arr(n_source)
      real*8 a_arr(n_source),b_arr(n_source)
      real*8 c_arr(n_source),d_arr(n_source)
      real*8 x_tmp(n_source+4),y_tmp(n_source+4)
      real*8 t(n_source+4),st(n_source+4)
      real*8 tmp1,tmp2
      integer i
c
c ------------------
c
      do i=1,n_source+4
         t(i) =0.d0
         st(i)=0.d0
      enddo
c
c  Calculate auxiliary points left and right.
c
      do i=1,n_source
         x_tmp(i)=x_arr(i)
         y_tmp(i)=y_arr(i)
      enddo
c
      CALL akima_ex(n_source,x_tmp,y_tmp,st)
c
c  Calculation of all polynomial slopes.
c
      do i=3,n_source+2
         if ( (st(i-2).eq.st(i-1)) .and. (st(i).eq.st(i+1)) ) then
            t(i) = 0.5d0 * (st(i+1)+st(i))
          else
            tmp1 = dABS( st(i+1) - st(i) )
            tmp2 = dABS( st(i-1) - st(i-2) )
            if ( (tmp1+tmp2).ne.0.d0 ) then
               t(i) = ( tmp1*st(i-1) + tmp2*st(i) ) / (tmp1+tmp2)
             else
               t(i) = 0.d0
            endif
         endif
      enddo
c
c  Calculation of polynomial coefficients.
c
      do i=1,n_source-1
c
         a_arr(i) = y_tmp(i+2)
         b_arr(i) = t(i+2)
c
         tmp1 = x_tmp(i+3) - x_tmp(i+2)
         if (tmp1.gt.1d-10) then
            c_arr(i) = ( ((3.0d0*st(i+2))-(2.d0*t(i+2)))-t(i+3) ) / tmp1
          else
            c_arr(i) = 0.d0
         endif
c
         tmp1 = tmp1*(x_tmp(i+3) - x_tmp(i+2))
         if (tmp1.gt.1d-10) THEN
            d_arr(i) = ( (t(i+2)+t(i+3)) - (2.d0*st(i+2)) ) / tmp1
           else
            d_arr(i) = 0.d0
         endif
c
      enddo
c
c ------------------
c
      return
      end
*
**************************************************************************
*
      real*8 function akima_po(n_source,x_arr,a,b,c,d,x)
c
!==========================================================================
!  Function   : akima_po
!  Input      : n_source  nr of function values
!               x_arr     abscissa values of tabulated function
!               a_arr, b_arr, c_arr, d_arr   constants of Akima interpolation
!               x         abscissa value for which function value is sought
!  Output     : interpolated/extrapolated function value at x
!  Description: Performes Akima interpolation at certain abscissa x
!  References : GOME level 1 extractor source code
!  Libraries  : none
!  Created    : 15-MAY-2001
!  Version nr : 1.0
!  Author     : Nick Schutgens, KNMI
!  Comments   : Originally by Akima. C implementation by Bernd Aberle,
!               DLR/WT-DA-SE from which this F90 code was developed
!--------------------------------------------------------------------------
c
      implicit none
c
c  Local variables
c
      integer n_source
      real*8 x_arr(n_source),x,delta
      real*8 a(n_source),b(n_source),c(n_source),d(n_source)
      integer klo,khi,k
c
c ------------------
c
      klo=1
      khi=n_source
      akima_po=0.d0
c
      do while (khi-klo.gt.1)
c
         k=(khi+klo)/2
         if (x_arr(k).gt.x) then
            khi=k
          else
            klo=k
         endif
         k=klo
c
         delta = x - x_arr(k)
         akima_po = a(k) + (b(k) * delta) + (c(k) * delta * delta) +
     #              (d(k) * delta * delta * delta)
c
      enddo
c
c ------------------
c
      return
      end
*
**************************************************************************