      SUBROUTINE COEINT(RAT,IR)
C ***************************************************************************
C
C
C            THE TRANSMISSION FUNCTION BETWEEN P1 AND P2 IS ASSUMED TO  HAVE
C       THE  FUNCTIONAL FORM
C                     TAU(P1,P2)= 1.0-SQRT(C*LOG(1.0+X*PATH)),
C               WHERE
C                     PATH(P1,P2)=((P1-P2)**2)*(P1+P2+CORE)/
C                                 (ETA*(P1+P2+CORE)+(P1-P2))
C
C
C        THE PARAMETERS C AND X ARE FUNCTIONS OF P2, AND ARE TO BE DETERMINED,
C        WHILE CORE IS A PRESPECIFIED NUMBER.ETA IS A FUNCTION OF THE THEIR
C        PRODUCT (CX);IT IS OBTAITED ITERATIVELY. THE DERIVATION OF ALL THESE
C        VALUES WILL BE EXPLAINED IN A FORTHCOMING PAPER.
C            SUBROUTINE COEINT DETERMINES C(I) AND X(I) BY USING THE ACTUAL
C        VALUES OF TAU(P(I-2),P(I)) AND TAU(P(I-1),P(I)) AND THE PREVIOUS
C        ITERATION VALUE OF ETA.
C             DEFINE: 
C                PATHA=PATH(P(I),P(I-2),CORE,ETA)
C                PATHB=PATH(P(I),P(I-1),CORE,ETA);
C        THEN
C                R=(1-TAU(P(I),P(I-2)))/(1-TAU(P(I),P(I-1)))
C                 = SQRT(LOG(1+X*PATHA)/LOG(1+X*PATHB)),
C        SO THAT
C                R**2= LOG(1+X*PATHA)/LOG(1+X*PATHB).
C        THIS EQUATION CAN BE SOLVED BY NEWTON'S METHOD FOR X AND THEN THE
C        RESULT USED TO FIND C. THIS IS REPEATED FOR EACH VALUE OF I GREATER
C        THAN 2 TO GIVE THE ARRAYS X(I) AND C(I).
C             NEWTON'S METHOD FOR SOLVING THE EQUATION
C                 F(X)=0
C        MAKES USE OF THE LOOP XNEW= XOLD-F(XOLD)/F'(XOLD).
C        THIS IS ITERATED 20 TIMES, WHICH IS PROBABLY EXCESSIVE.
C        THE FIRST GUESS FOR ETA IS 3.2E-4*EXP(-P(I)/1000),WHICH HAS
C        BEEN FOUND TO BE FAIRLY REALISTIC BY EXPERIMENT; WE ITERATE 5 TIMES
C        (AGAIN,PROBABLY EXCESSIVELY) TO OBTAIN THE VALUES FOR C,X,ETA TO BE
C        USED FOR INTERPOLATION.
C           THERE ARE SEVERAL POSSIBLE PITFALLS: 
C              1) IN THE COURSE OF ITERATION, X MAY REACH A VALUE WHICH MAKES
C                 1+X*PATHA NEGATIVE; IN THIS CASE THE ITERATION IS STOPPED,
C                 AND AN ERROR MESSAGE IS PRINTED OUT.
C              2) EVEN IF (1) DOES NOT OCCUR, IT IS STILL POSSIBLE THAT X MAY
C                 BE NEGATIVE AND LARGE ENOUGH TO MAKE 1+X*PATH(P(I),0,CORE,ETA)
C                 NEGATIVE. THIS IS CHECKED FOR IN A FINAL LOOP, AND IF TRUE,
C                 A WARNING IS PRINTED OUT.
C
C  **************************************************************************
      implicit real(selected_real_kind(15)) (a-h,o-z)
      COMMON/PRESS/PA(109)
      REAL(selected_real_kind(15)) :: RAT
      REAL(selected_real_kind(15)) :: PA,TRANSA
      REAL(selected_real_kind(15)) :: PA2
      COMMON/TRAN/ TRANSA(109,109)
      COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP
      DIMENSION PATH0(109),ETAP(109),XAP(109),CAP(109)
      DIMENSION SINV(4)
      DATA SINV/2.74992,2.12731,4.38111,.0832926/
      CORE=5.000
      UEXP=0.90
      P0=0.7
      DO 902 I=1,109
      PA2=PA(I)*PA(I)
      SEXPV(I)=.505+2.0E-5*PA(I)+.035*(PA2-.25)/(PA2+.25)
902   CONTINUE
      DO 900 I=1,109
      ETA(I)=3.2E-4*EXP(-PA(I)/500.)
      ETAP(I)=ETA(I)
900   CONTINUE
      DO 1200 NP=1,10
      DO 1000 I=3,109
      padi=pa(i)
      padim=pa(i-1)
      padim2=pa(i-2)
      SEXP=SEXPV(I)
      R=(1.0D0-TRANSA(I,I-2))/(1.0D0-TRANSA(I,I-1))
      REXP=R**(UEXP/SEXP)
      PATHA=(PATH(padi,padim2,CORE,ETA(I)))**UEXP
      PATHB=(PATH(padi,padim,CORE,ETA(I)))**UEXP
      XX=2.0D0*(PATHB*REXP-PATHA)/(PATHB*PATHB*REXP-PATHA*PATHA)
      DO 1010 LL=1,20
      ftest1=xx*patha
      ftest2=xx*pathb
c*** end iteration and solve if ftest1 is small or ftest2 is large
      if (ftest1.le.1.0d-10) then
          xx=1.0d0
          xa(i)=xx
          ca(i)=(1.0d0-transa(i,i-2))**(uexp/sexp)/patha
          go to 1011
      endif
      if (ftest2.ge.1.0d8) then
          xxlog=(dlog(patha)-rexp*dlog(pathb))/(rexp-1.0d0)
          xx=dlog(xxlog)
          xa(i)=xx
          ca(i)=(1.0d0-transa(i,i-2))**(uexp/sexp)/(xxlog+dlog(patha))
          go to 1011
      endif
      F1=DLOG(1.0D0+XX*PATHA)
      F2=DLOG(1.0D0+XX*PATHB)
      F=F1/F2-REXP
      FPRIME=(F2*PATHA/(1.0D0+XX*PATHA)-F1*PATHB/(1.0D0+XX*PATHB))/
     1    (F2*F2)
      XX=XX-F/FPRIME
c     Added this fixup to ensure convergence mrd 8/12/92
      if ( xx*patha .le. -1 ) then
	  print*, ' Fixup ', xx, -(1.-1.d-10)/patha
	  xx = -(1.-1.d-10)/patha
      end if
      CHECK=1.0D0+XX*PATHA
      IF (CHECK) 1020,1020,1025
 1020 CONTINUE
      WRITE (6,360) I,LL,CHECK
360   FORMAT (' ERROR,I=',I3,'LL=',I3,'CHECK=',F20.10)
      STOP
1025  CONTINUE
 1010 CONTINUE
      CA(I)=(1.0D0-TRANSA(I,I-2))**(UEXP/SEXP)/
     1 (DLOG(1.0D0+XX*PATHA)+1.0D-20)
      XA(I)=XX
 1011  continue
1000  CONTINUE
      XA(2)=XA(3)
      XA(1)=XA(3)
      CA(2)=CA(3)
      CA(1)=CA(3)
      DO 1100 I=3,109
      padi=pa(i)
!     This change is necessary for atmos.
      PATH0(I)=(PATH(padi,0.0d0,CORE,ETA(I)))**UEXP
!      PATH0(I)=(PATH(padi,0.,CORE,ETA(I)))**UEXP
      PATH0(I)=1.0D0+XA(I)*PATH0(I)
      IF (PATH0(I).LT.0.) WRITE (6,361) I,PATH0(I),XA(I)
1100  CONTINUE
      DO 1035 I=1,109
      drat=rat
      SEXP=SEXPV(I)
      ETAP(I)=ETA(I)
      ETA(I)=(SINV(IR)/drat)**(1./SEXP)*
     1  (CA(I)*XA(I))**(1./UEXP)
1035  CONTINUE
C
C     THE ETA FORMULATION IS DETAILED IN SCHWARZKOPF AND FELS(1985).
C        THE QUANTITY SINV=(G*DELTANU)/(RCO2*D*S)
C      IN CGS UNITS,WITH D,THE DIFFUSICITY FACTOR=2, AND
C      S,THE SUM OF CO2 LINE STRENGTHS OVER THE 15UM CO2 BAND
C       ALSO,THE DENOMINATOR IS MULTIPLIED BY
C      1000 TO PERMIT USE OF MB UNITS FOR PRESSURE.
C        S IS ACTUALLY WEIGHTED BY B(250) AT 10 CM-1 WIDE INTERVALS,IN
C      ORDER TO BE CONSISTENT WITH THE METHODS USED TO OBTAIN THE LBL
C      1-BAND CONSOLIDATED TRANCMISSION FUNCTIONS.
C      FOR THE 490-850 INTERVAL (DELTANU=360,IR=1) SINV=2.74992.
C      (SLIGHTLY DIFFERENT FROM 2.7528 USED IN EARLIER VERSIONS)
C      FOR THE 490-670 INTERVAL (IR=2) SINV=2.12731
C      FOR THE 670-850 INTERVAL (IR=3) SINV=4.38111
C      FOR THE 2270-2380 INTERVAL (IR=4) SINV=0.0832926
C      SINV HAS BEEN OBTAINED USING THE 1982 AFGL CATALOG FOR CO2
C        RAT IS THE ACTUAL CO2 MIXING RATIO IN UNITS OF 330 PPMV,
C      LETTING USE OF THIS FORMULATION FOR ANY CO2 CONCENTRATION.
C
C     WRITE (6,366) (NP,I,CA(I),XA(I),ETA(I),SEXPV(I),I=1,109)
C366   FORMAT (2I4,4E20.12)
1200  CONTINUE
 361  FORMAT (' **WARNING:** 1+XA*PATH(PA(I),0) IS NEGATIVE,I= ',I3,/
     1 20X,'PATH0(I)=',F16.6,' XA(I)=',F16.6)
      RETURN
      END
