#include "dims.h"
      SUBROUTINE RECUR
      implicit none
#include "params.h"
#include "radcool.h"
!
! Local:
      real :: CO2INT(4),UREF(4),A(zimxp),COR(zimxp),UC(zimxp)
!
      real :: amat,bmat,al
      COMMON /CO2CFG/ AMAT(43,9),BMAT(43,9)
      common /CO2CFG_PVT/ AL(zimxp,nzpm9_zpm5)
!$OMP THREADPRIVATE (/co2cfg_pvt/)
!DIR$ TASKCOMMON CO2CFG_PVT
      real :: xr,xrh,xrf
      COMMON/PIRGRD/ XR(nzpsrf_zpm5),XRH(59),XRF(nzpm9_zpm5)
      real :: uco2ro,alo,cor150,cor360,cor540,cor720,uco2co
      common /PIRNE/ uco2ro(51),alo(51),
     1               cor150(6),cor360(6),cor540(6),cor720(6),uco2co(6)
!
! Local:
      integer :: len1,k,i
C
C     ****  UCO2(17) (CO2 COLUMN AMOUNT) FOR CO2
C
C     **** CALCULATE COEFICIENTS FOR THE RECCURENCE FORMULA:
C
C     **** BETWEEN X=12.5 AND 13.75 THESE COEFFICIENTS (AL) ARE
C     **** CALCULATED USING CORRECTIONS TO THE ESCAPE FUNCTION.
C     **** STARTING FROM X=14.00 AND ABOVE THE PARAMETERIZATION
C     **** COEFFICIENTS ARE EQUAL TO THE ESCAPE FUNCTION.
C
      LEN1 = zimxp
      DO 3 K=1,6
        CO2INT(1) = COR150(K)
        CO2INT(2) = COR360(K)
        CO2INT(3) = COR540(K)
        CO2INT(4) = COR720(K)
        UREF(1) = UCO2CO(K)*150./360.
        UREF(2) = UCO2CO(K)
        UREF(3) = UCO2CO(K)*540./360.
        UREF(4) = UCO2CO(K)*720./360.
        DO 4 I=1,LEN1
          UC(I) = UCO2(I,K)
  4     CONTINUE
        CALL A18LINV(UC,A,UCO2RO,ALO,51,LEN1)
        CALL A18LINV(UC,COR,UREF,CO2INT,4,LEN1)
        DO 5 I=1,LEN1
          AL(I,K) = EXP(COR(I)+A(I))
  5     CONTINUE
  3   CONTINUE
      DO 6 K=7,nzpm9_zpm5
        DO 7 I=1,LEN1
          UC(I) = UCO2(I,K)
  7     CONTINUE
        CALL A18LINV(UC,A,UCO2RO,ALO,51,LEN1)
        DO 8 I=1,LEN1
          AL(I,K) = EXP(A(I))
  8     CONTINUE
  6   CONTINUE
      RETURN
      END
C
      SUBROUTINE A18LINV(X,Y,XN,YN,N,IMAX)
      implicit none
C     ****
C     ****     This procedure performs linear interpolation within the
C     ****     table defined by the N points (XN(NN),Y(NN)).
C     ****     Where:
C     ****
C     ****       NN = 1,N,1
C     ****
C     ****       XN(NN) < XN(NN+1) for NN = 1,N-1
C     ****
C     ****     Parameters:
C     ****
C     ****       X(IMAX) = array of IMAX x-values at which linear
C     ****                 interpolation is required
C     ****
C     ****       XN(N) = array of N abscissae at which function values
C     ****               are given
C     ****
C     ****       YN(N) = function values corresponding to abscissae,
C     ****               XN(N)
C     ****
C     ****     Output:
C     ****
C     ****      Y(IMAX)  The IMAX interpolated values are
C     ****                     returned in this array
C     ****
!
! Args:
      integer,intent(in) :: n,imax
      real,intent(out)   :: Y(IMAX)
      real,intent(in)    :: X(IMAX), XN(N), YN(N) 
!
! Local:
      integer :: KK(IMAX),i,nn
C     ****
C     ****     Where:
C     ****       Y(IMAX) is vector output
C     ****
C     ****       KK is work space
C     ****
C     ****
C     ****     Initialize array KK
C     ****
      DO I = 1,IMAX
	KK(I) = 0
      ENDDO
C     ****
C     ****     Locate interval in (XN,YN) in table containing X(I)
C     ****
      DO NN = 1,N-1
	DO I = 1,IMAX
!         KK(I) = merge(NN+1,KK(I),(XN(NN+1)-X(I))*(X(I)-XN(NN))>=0.)
          if ((xn(nn+1)-x(i))*(x(i)-xn(nn)) >= 0.) kk(i) = nn+1
	ENDDO
      ENDDO
C     ****
C     ****     Check for
C     ****
C     ****       X(I) < XN(1),  X(I) > X(N)
C     ****
C     ****       and use linear extrapolation if necessary
C     ****
      DO I = 1,IMAX
!       KK(I) = merge(2,KK(I),XN(1)-X(I)>=0.)
!       KK(I) = merge(N,KK(I),X(I)-XN(N)>=0.)
        if (xn(1)-x(i) >= 0.) kk(i) = 2
        if (x(i) >= xn(n)) kk(i) = n
      ENDDO
C     ****
C     ****     Perform interpolation prescribed above
C     ****
      DO I = 1,IMAX
	Y(I) = (YN(KK(I)-1)*(XN(KK(I))-X(I)) + YN(KK(I))*
     1               (X(I)-XN(KK(I)-1)))/(XN(KK(I))-XN(KK(I)-1))
      ENDDO
      RETURN
      END
C