      FUNCTION ANTEMP(L,Z)
      REAL :: ZB(10,7),C(11,7),DELTA(10,7),TSTAR(7)
! ************** TROPICAL SOUNDING **************************
      DATA (ZB(N,1),N=1,10)/  2.0,   3.0,   16.5,  21.5,  45.0, &
                              51.0,  70.0,  100.,  200.,  300./
      DATA (C(N,1),N=1,11)/ -6.0,  -4.0,  -6.7,   4.0,   2.2,   &
                         1.0,  -2.8,  -.27,   0.0,   0.0,  0.0/
      DATA (DELTA(N,1),N=1,10)/.5,    .5,    .3,    .5,    1.0, &
                              1.0,   1.0,   1.0,   1.0,    1.0/
! ************** SUB-TROPICAL SUMMER ************************
      DATA (ZB(N,2),N=1,10)/ 1.5,   6.5,  13.0,  18.0,  26.0, &
                              36.0,  48.0,  50.0, 70.0,  100./
      DATA (C(N,2),N=1,11)/ -4.0,  -6.0,  -6.5,   0.0,   1.2, &
                        2.2,   2.5,   0.0,  -3.0,  -0.25,  0.0/
      DATA (DELTA(N,2),N=1,10)/ .5,  1.0,    .5,    .5,   1.0, &
                              1.0,  2.5,    .5,   1.0,   1.0/
! ************** SUB-TROPICAL WINTER ************************
      DATA (ZB(N,3),N=1,10)/ 3.0,  10.0,  19.0,  25.0,  32.0, &
                              44.5, 50.0,  71.0,  98.0,  200.0/
      DATA (C(N,3),N=1,11)/ -3.5,  -6.0,  -0.5,  0.0,   0.4, &
                              3.2,   1.6,  -1.8, -0.7,   0.0,   0.0/
      DATA (DELTA(N,3),N=1,10)/ .5,   .5,  1.0,   1.0,   1.0, &
                              1.0,  1.0,  1.0,   1.0,   1.0/
! *************  SUB-ARCTIC SUMMER *************************
      DATA (ZB(N,4),N=1,10)/ 4.7, 10.0,  23.0,  31.8,  44.0, &
                              50.2, 69.2, 100.0, 102.0, 103.0/
      DATA (C(N,4),N=1,11)/ -5.3, -7.0,   0.0,  1.4,   3.0, &
                               0.7, -3.3,  -0.2,  0.0,   0.0,  0.0/
      DATA (DELTA(N,4),N=1,10)/ .5,   .3,  1.0,   1.0,   2.0, &
                              1.0,  1.5,  1.0,   1.0,   1.0/
! ************ SUB-ARCTIC WINTER *****************************
      DATA (ZB(N,5),N=1,10)/ 1.0,   3.2,   8.5,   15.5,   25.0, &
                              30.0,  35.0,  50.0,  70.0,  100.0/
      DATA (C(N,5),N=1,11)/ 3.0,  -3.2,  -6.8,  0.0,  -0.6, &
                              1.0,   1.2,   2.5, -0.7,  -1.2,  0.0/
      DATA (DELTA(N,5),N=1,10)/ .4,   1.5,    .3 ,   .5,   1.0, &
                              1.0,   1.0,   1.0,   1.0,   1.0/
! ************ US STANDARD 1976 ******************************
      DATA (ZB(N,6),N=1,10)/ 11.0,  20.0,  32.0,  47.0,  51.0, & 
                             71.0,  84.8520,  90.0,  91.0,  92.0/
      DATA (C(N,6),N=1,11)/ -6.5,   0.0,   1.0,   2.80,  0.0, &
                             -2.80, -2.00,  0.0,   0.0,   0.0,  0.0/
      DATA (DELTA(N,6),N=1,10)/ 0.3,   1.0,   1.0,   1.0,   1.0, &
                              1.0,   1.0,   1.0,   1.0,   1.0/
!
! ************ ENLARGED US STANDARD 1976 **********************
      DATA (ZB(N,7),N=1,10)/ 11.0,  20.0,  32.0,  47.0,  51.0, &
                             71.0,  84.8520,  90.0,  91.0,  92.0/
      DATA (C(N,7),N=1,11)/ -6.5,   0.0,   1.0,   2.80,  0.0, &
                             -2.80, -2.00,  0.0,   0.0,   0.0,  0.0/
      DATA (DELTA(N,7),N=1,10)/ 0.3,   1.0,   1.0,   1.0,   1.0, &
                              1.0,   1.0,   1.0,   1.0,   1.0/
!
      DATA TSTAR/ 300.0,  294.0,  272.2,  287.0,  257.1, 2*288.15/
!
      NLAST=10
      TEMP=TSTAR(L)+C(1,L)*Z
      DO 20 N=1,NLAST
      EXPO=(Z-ZB(N,L))/DELTA(N,L)
      EXPP=ZB(N,L)/DELTA(N,L)
!JD single-precision change
!      FAC=EXP(EXPP)+EXP(-EXPP)
!mp	write(6,*) '.........................................'
!mp what in the hell does the next line do?
!mp	
!mp	apparently if statement <0 or =0 then 23, else 24
!mp     IF(ABS(EXPO)-100.0) 23,23,24
!
! changed to a more reasonable value for the workstation	
!
      IF(ABS(EXPO)-50.0) 23,23,24
   23 X=EXP(EXPO)
      Y=X+1.0/X
      ZLOG=ALOG(Y)
      GO TO 25
   24 ZLOG=ABS(EXPO)
!mp   25 IF(EXPP-100.0) 27,27,28
   25 IF(EXPP-50.0) 27,27,28
!JD single-precision change
   27 FAC=EXP(EXPP)+EXP(-EXPP)
      FACLOG=ALOG(FAC)
      GO TO 29
   28 FACLOG=EXPP
!     TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)*
!    1     ALOG((EXP(EXPO)+EXP(-EXPO))/FAC))
   29 TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)* &
           (ZLOG-FACLOG))
!mp	write(6,*) 'ANTEMP pieces (C,C,ZLOG,FACLOG)', C(N+1,L),C(N,L),
!mp     +	ZLOG,FACLOG
   20 CONTINUE
      ANTEMP=TEMP

      END FUNCTION ANTEMP

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc

      SUBROUTINE COEINT(RAT,IR)
! **********************************************************************
!
!
!            THE TRANSMISSION FUNCTION BETWEEN P1 AND P2 IS ASSUMED TO
!       THE  FUNCTIONAL FORM
!                     TAU(P1,P2)= 1.0-SQRT(C*LOG(1.0+X*PATH)),
!               WHERE
!                     PATH(P1,P2)=((P1-P2)**2)*(P1+P2+CORE)/
!                                 (ETA*(P1+P2+CORE)+(P1-P2))
!
!
!        THE PARAMETERS C AND X ARE FUNCTIONS OF P2, AND ARE TO BE DETER
!        WHILE CORE IS A PRESPECIFIED NUMBER.ETA IS A FUNCTION OF THE TH
!        PRODUCT (CX);IT IS OBTAITED ITERATIVELY. THE DERIVATION OF ALL
!        VALUES WILL BE EXPLAINED IN A FORTHCOMING PAPER.
!            SUBROUTINE COEINT DETERMINES C(I) AND X(I) BY USING THE ACT
!        VALUES OF TAU(P(I-2),P(I)) AND TAU(P(I-1),P(I)) AND THE PREVIOU
!        ITERATION VALUE OF ETA.
!             DEFINE:
!                PATHA=PATH(P(I),P(I-2),CORE,ETA)
!                PATHB=PATH(P(I),P(I-1),CORE,ETA);
!        THEN
!                R=(1-TAU(P(I),P(I-2)))/(1-TAU(P(I),P(I-1)))
!                 = SQRT(LOG(1+X*PATHA)/LOG(1+X*PATHB)),
!        SO THAT
!                R**2= LOG(1+X*PATHA)/LOG(1+X*PATHB).
!        THIS EQUATION CAN BE SOLVED BY NEWTON S METHOD FOR X AND THEN T
!        RESULT USED TO FIND C. THIS IS REPEATED FOR EACH VALUE OF I GRE
!        THAN 2 TO GIVE THE ARRAYS X(I) AND C(I).
!             NEWTON S METHOD FOR SOLVING THE EQUATION
!                 F(X)=0
!        MAKES USE OF THE LOOP XNEW= XOLD-F(XOLD)/FPRIME(XOLD).
!        THIS IS ITERATED 20 TIMES, WHICH IS PROBABLY EXCESSIVE.
!        THE FIRST GUESS FOR ETA IS 3.2E-4*EXP(-P(I)/1000),WHICH HAS
!        BEEN FOUND TO BE FAIRLY REALISTIC BY EXPERIMENT; WE ITERATE 5 T
!        (AGAIN,PROBABLY EXCESSIVELY) TO OBTAIN THE VALUES FOR C,X,ETA T
!        USED FOR INTERPOLATION.
!           THERE ARE SEVERAL POSSIBLE PITFALLS:
!              1) IN THE COURSE OF ITERATION, X MAY REACH A VALUE WHICH
!                 1+X*PATHA NEGATIVE; IN THIS CASE THE ITERATION IS STOP
!                 AND AN ERROR MESSAGE IS PRINTED OUT.
!              2) EVEN IF (1) DOES NOT OCCUR, IT IS STILL POSSIBLE THAT
!                 BE NEGATIVE AND LARGE ENOUGH TO MAKE 1+X*PATH(P(I),0,C
!                 NEGATIVE. THIS IS CHECKED FOR IN A FINAL LOOP, AND IF
!                 A WARNING IS PRINTED OUT.
!
!  *********************************************************************
!....
!     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!     COMMON/PRESS/PA(109)
      REAL RAT,SINV
!     REAL PA,CORE,TRANSA,PATH,UEXP,SEXP,ETA,SEXPV
      REAL 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)
      INTEGER :: IERR
      DATA SINV/2.74992,2.12731,4.38111,0.0832926/
!NOV89   DIMENSION SINV(3)
!NOV89   DATA SINV/2.74992,2.12731,4.38111/
!O222  OLD CODE USED 2.7528 RATHER THAN 2.74992 ---K.A.C. OCTOBER 1988
!O222   WHEN 2.7528 WAS USED,WE EXACTLY REPRODUCED THE MRF CO2 ARRAYS
      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
      SEXP=SEXPV(I)
      R=(1.0D0-TRANSA(I,I-2))/(1.0D0-TRANSA(I,I-1))
      REXP=R**(UEXP/SEXP)
      arg1=path(pa(i),pa(i-2),core,eta(i))
      arg2=path(pa(i),pa(i-1),core,eta(i))
      PATHA=(PATH(PA(I),PA(I-2),CORE,ETA(I)))**UEXP
      PATHB=(PATH(PA(I),PA(I-1),CORE,ETA(I)))**UEXP
      XX=2.0D0*(PATHB*REXP-PATHA)/(PATHB*PATHB*REXP-PATHA*PATHA)
      DO 1010 LL=1,20
      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))/ &
          (F2*F2)
      XX=XX-F/FPRIME
      CHECK=1.0D0+XX*PATHA
!!!!  IF (CHECK) 1020,1020,1025
      IF(CHECK.LE.0.)THEN
        WRITE(errmess,360)I,LL,CHECK
        WRITE(errmess,*)' xx=',xx,' patha=',patha
  360   FORMAT(' ERROR,I=',I3,'LL=',I3,'CHECK=',F20.10)
        CALL wrf_error_fatal ( errmess )
      ENDIF
 1010 CONTINUE
      CA(I)=(1.0D0-TRANSA(I,I-2))**(UEXP/SEXP)/ &
       (DLOG(1.0D0+XX*PATHA)+1.0D-20)
      XA(I)=XX
1000  CONTINUE
      XA(2)=XA(3)
      XA(1)=XA(3)
      CA(2)=CA(3)
      CA(1)=CA(3)
      DO 1100 I=3,109
      PATH0(I)=(PATH(PA(I),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
      SEXP=SEXPV(I)
      ETAP(I)=ETA(I)
      ETA(I)=(SINV(IR)/RAT)**(1./SEXP)* &
        (CA(I)*XA(I))**(1./UEXP)
1035  CONTINUE
!
!     THE ETA FORMULATION IS DETAILED IN SCHWARZKOPF AND FELS(1985).
!        THE QUANTITY SINV=(G*DELTANU)/(RCO2*D*S)
!      IN CGS UNITS,WITH D,THE DIFFUSICITY FACTOR=2, AND
!      S,THE SUM OF CO2 LINE STRENGTHS OVER THE 15UM CO2 BAND
!       ALSO,THE DENOMINATOR IS MULTIPLIED BY
!      1000 TO PERMIT USE OF MB UNITS FOR PRESSURE.
!        S IS ACTUALLY WEIGHTED BY B(250) AT 10 CM-1 WIDE INTERVALS,IN
!      ORDER TO BE CONSISTENT WITH THE METHODS USED TO OBTAIN THE LBL
!      1-BAND CONSOLIDATED TRANCMISSION FUNCTIONS.
!      FOR THE 490-850 INTERVAL (DELTANU=360,IR=1) SINV=2.74992.
!      (SLIGHTLY DIFFERENT FROM 2.7528 USED IN EARLIER VERSIONS)
!      FOR THE 490-670 INTERVAL (IR=2) SINV=2.12731
!      FOR THE 670-850 INTERVAL (IR=3) SINV=4.38111
!      FOR THE 2270-2380 INTERVAL (IR=4) SINV=0.0832926
!      SINV HAS BEEN OBTAINED USING THE 1982 AFGL CATALOG FOR CO2
!        RAT IS THE ACTUAL CO2 MIXING RATIO IN UNITS OF 330 PPMV,
!      LETTING USE OF THIS FORMULATION FOR ANY CO2 CONCENTRATION.
!
!     WRITE (6,366) (NP,I,CA(I),XA(I),ETA(I),SEXPV(I),I=1,109)
!366   FORMAT (2I4,4E20.12)
1200  CONTINUE
 361  FORMAT (' **WARNING:** 1+XA*PATH(PA(I),0) IS NEGATIVE,I= ',I3,/ &
       20X,'PATH0(I)=',F16.6,' XA(I)=',F16.6)
      RETURN
      END SUBROUTINE COEINT

!--------------


!CCC  PROGRAM CO2INS
      SUBROUTINE CO2INS(T22,T23,T66,IQ,L,LP1,iflag)
!     *********************************************************
!       SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 ******
!          ..... K.CAMPANA   MARCH 1988,OCTOBER 1988...
!          ..... K.CAMPANA   DECEMBER 1988-CLEANED UP FOR LAUNCHER
!          ..... K.CAMPANA   NOVEMBER 1989-ALTERED FOR NEW RADIATION
!     *********************************************************
      DIMENSION T22(LP1,LP1,3),T23(LP1,LP1,3),T66(LP1,LP1,6)
      DIMENSION DCDT8(LP1,LP1),DCDT10(LP1,LP1),CO2PO(LP1,LP1), &
       CO2800(LP1,LP1),CO2PO1(LP1,LP1),CO2801(LP1,LP1),CO2PO2(LP1,LP1), &
       CO2802(LP1,LP1),N(LP1),D2CT8(LP1,LP1),D2CT10(LP1,LP1)
!CC   ITIN=22
!CC   ITIN1=23
!O222  LATEST CODE HAD  IQ=1
!CC      IQ=4
1011  FORMAT (4F20.14)
!CC      READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1)
!CC      READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1)
!CC      READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1)
!CC      READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1)
!CC      READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1)
!CC      READ (ITIN1,1011) ((CO2802(I,J),I=1,LP1),J=1,LP1)
      DO 300 J=1,LP1
        DO 300 I=1,LP1
          CO2PO(I,J) = T22(I,J,1)
!NOV89
          IF (IQ.EQ.5) GO TO 300
!NOV89
          CO2PO1(I,J) = T22(I,J,2)
          CO2PO2(I,J) = T22(I,J,3)
  300 CONTINUE
      DO 301 J=1,LP1
        DO 301 I=1,LP1
          CO2800(I,J) = T23(I,J,1)
!NOV89
          IF (IQ.EQ.5) GO TO 301
!NOV89
          CO2801(I,J) = T23(I,J,2)
          CO2802(I,J) = T23(I,J,3)
  301 CONTINUE
!***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS
!   ARE:
!        IQ=1    560-800     (CONSOL.=490-850)
!        IQ=2    560-670     (CONSOL.=490-670)
!        IQ=3    670-800     (CONSOL.=670-850)
!        IQ=4    560-760 (ORIGINAL CODE)   (CONSOL.=490-850)
!NOV89
!        IQ=5   2270-2380    (CONSOL.=2270-2380)
!NOV89
!  THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS
!  USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT
!  WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S.
!NOV89
!      NOTE: ALTHOUGH THE BAND TRANSMISSION FUNCTIONS ARE
!  COMPUTED FOR ALL RADIATIVE BANDS, AS OF 9/28/88, THEY
!  ARE WRITTEN OUT IN FULL ONLY FOR THE FULL 15 UM BAND CASES
!  (IQ=1,4).  IN OTHER CASES, THE TRANSMISSIVITIES (1,K) ARE
!  WRITTEN OUT, AS THESE ARE THE ONLY ONES NEEDED FOR CTS
!  CALCULATIONS.  ALSO, FOR THE 4.3 UM BAND (IQ=5) THE TEMP.
!  DERIVATIVE TERMS ARE NOT WRITTEN OUT, AS THEY ARE UNUSED.
!NOV89
      IF (IQ.EQ.1) THEN
         C1=1.5
         C2x=0.5
      ENDIF
      IF (IQ.EQ.2) THEN
        C1=18./11.
        C2x=7./11.
      ENDIF
      IF (IQ.EQ.3) THEN
        C1=18./13.
        C2x=5./13.
      ENDIF
      IF (IQ.EQ.4) THEN
        C1=1.8
        C2x=0.8
      ENDIF
!NOV89
      IF (IQ.EQ.5) THEN
        C1=1.0
        C2x=0.0
      ENDIF
!NOV89
      DO 1021 I=1,LP1
      DO 1021 J=1,LP1
      CO2PO(J,I)=C1*CO2PO(J,I)-C2x
      CO2800(J,I)=C1*CO2800(J,I)-C2x
!NOV89
      IF (IQ.EQ.5) GO TO 1021
!NOV89
      CO2PO1(J,I)=C1*CO2PO1(J,I)-C2x
      CO2801(J,I)=C1*CO2801(J,I)-C2x
      CO2PO2(J,I)=C1*CO2PO2(J,I)-C2x
      CO2802(J,I)=C1*CO2802(J,I)-C2x
1021  CONTINUE
!NOV89
      IF (IQ.GE.1.AND.IQ.LE.4) THEN
!NOV89
      DO 1 J=1,LP1
      DO 1 I=1,LP1
      DCDT8(I,J)=.02*(CO2801(I,J)-CO2802(I,J))*100.
      DCDT10(I,J)=.02*(CO2PO1(I,J)-CO2PO2(I,J))*100.
      D2CT8(I,J)=.0016*(CO2801(I,J)+CO2802(I,J)-2.*CO2800(I,J))*1000.
      D2CT10(I,J)=.0016*(CO2PO1(I,J)+CO2PO2(I,J)-2.*CO2PO(I,J))*1000.
1     CONTINUE
!NOV89
      ENDIF
!NOV89
!O222 *********************************************************
!CC       REWIND 66
!        SAVE CDT51,CO251,C2D51,CDT58,CO258,C2D58..ON TEMPO FILE
!CC       WRITE (66) DCDT10
!CC       WRITE (66) CO2PO
!CC       WRITE (66) D2CT10
!CC       WRITE (66) DCDT8
!CC       WRITE (66) CO2800
!CC       WRITE (66) D2CT8
!CC       REWIND 66
!NOV89
      IF (IQ.EQ.1.OR.IQ.EQ.4) THEN
!NOV89
      DO 400 J=1,LP1
       DO 400 I=1,LP1
        T66(I,J,1) = DCDT10(I,J)
        T66(I,J,2) = CO2PO(I,J)
        T66(I,J,3) = D2CT10(I,J)
        T66(I,J,4) = DCDT8(I,J)
        T66(I,J,5) = CO2800(I,J)
        T66(I,J,6) = D2CT8(I,J)
  400 CONTINUE
!NOV89
      ELSE
      DO 409 I=1,LP1
        T66(I,1,2) = CO2PO(1,I)
        T66(I,1,5) = CO2800(1,I)
        IF (IQ.EQ.5) GO TO 409
        T66(I,1,1) = DCDT10(1,I)
        T66(I,1,3) = D2CT10(1,I)
        T66(I,1,4) = DCDT8(1,I)
        T66(I,1,6) = D2CT8(1,I)
  409 CONTINUE
      ENDIF
!NOV89
!O222 *********************************************************
      RETURN
      END SUBROUTINE CO2INS
!O222 PROGRAM CO2INT(INPUT,TAPE5=INPUT)
!NOV89
      SUBROUTINE CO2INT(ITAPE,T15A,T15B,T22,RATIO,IR,NMETHD,NLEVLS,NLP1,NLP2)
!NOV89
!     *********************************************************
!       CHANGES TO DATA READ  AND FORMAT SEE CO222     ***
!          ..... K.CAMPANA   MARCH 1988,OCTOBER 1988
!       CHANGES TO PASS ITAPE,AND IF IR=4,READ 1 CO2 REC..KAC NOV89
!     *********************************************************
!       CO2INT INTERPOLATES CARBON DIOXIDE TRANSMISSION FUNCTIONS
!  FROM THE 109 LEVEL GRID,FOR WHICH THE TRANSMISSION FUNCTIONS
!  HAVE BEEN PRE-CALCULATED, TO THE GRID STRUCTURE SPECIFIED BY THE
!  USER.
!
!        METHOD:
!
!      CO2INT IS EMPLOYABLE FOR TWO PURPOSES: 1) TO OBTAIN TRANSMIS-
!  SIVITIES BETWEEN ANY 2 OF AN ARRAY OF USER-DEFINED PRESSURES; AND
!  2) TO OBTAIN LAYER-MEAN TRANSMISSIVITIES BETWEEN ANY 2 OF AN ARRAY
!  OF USER-DEFINED PRESSURE LAYERS.TO CLARIFY THESE TWO PURPOSES,SEE
!  THE DIAGRAM AND DISCUSSION BELOW.
!      CO2INT MAY BE USED TO EXECUTE ONLY ONE PURPOSE AT ONE TIME.
!
!     LET P BE AN ARRAY OF USER-DEFINED PRESSURES
!     AND PD BE USER-DEFINED PRESSURE LAYERS.
!
!       - - - - - - - - -   PD(I-1) ---
!                                     ^
!       -----------------   P(I)      ^  PRESSURE LAYER I  (PLM(I))
!                                     ^
!       - - - - - - - - -   PD(I)  ---
!                                     ^
!       -----------------   P(I+1)    ^  PRESSURE LAYER I+1 (PLM(I+1))
!                                     ^
!       - - - - - - - - -   PD(I+1)---
!            ...                          (THE NOTATION USED IS
!            ...                          CONSISTENT WITH THE CODE)
!            ...
!      - - - - - - - - -    PD(J-1)
!
!      -----------------    P(J)
!
!      - - - - - - - - -    PD(J)
!
!      PURPOSE 1:   THE TRANSMISSIVITY BETWEEN SPECIFIC PRESSURES
!      P(I) AND P(J) ,TAU(P(I),P(J))  IS COMPUTED BY THIS PROGRAM.
!      IN THIS MODE,THERE IS NO REFERENCE TO LAYER PRESSURES PD
!      (PD,PLM ARE NOT INPUTTED).
!
!      PURPOSE 2:   THE LAYER-MEAN TRANSMISSIVITY BETWEEN A LAYER-
!      MEAN PRESSURE PLM(J) AND PRESSURE LAYER I IS GIVEN BY
!         TAULM(PLM(I),PLM(J)). IT IS COMPUTED BY THE INTEGRAL
!
!                           PD(I)
!                           ----
!             1             ^
!        -------------  *   ^   TAU ( P',PLM(J) )  DP'
!        PD(I)-PD(I-1)      ^
!                        ----
!                        PD(I-1)
!
!           THE LAYER-MEAN PRESSURE PLM(I) IS SPECIFIED BY THE USER.
!        FOR MANY PURPOSES,PLM WILL BE CHOSEN TO BE THE AVERAGE
!        PRESSURE IN THE LAYER-IE,PLM(I)=0.5*(PD(I-1)+PD(I)).
!           FOR LAYER-MEAN TRANSMISSIVITIES,THE USER THUS INPUTS
!        A PRESSURE ARRAY (PD) DEFINING THE PRESSURE LAYERS AND AN
!        ARRAY (PLM) DEFINING THE LAYER-MEAN PRESSURES.THE CALCULATION
!        DOES NOT DEPEND ON THE P ARRAY USED FOR PURPOSE 1 (P IS NOT
!        INPUTTED).
!
!            THE FOLLOWING PARAGRAPHS DEPICT THE UTILIZATION OF THIS
!       CODE WHEN USED TO COMPUTE TRANSMISSIVITIES BETWEEN SPECIFIC
!       PRESSURES. LATER PARAGRAPHS DESCRIBE ADDITIONAL FEATURES NEEDED
!       FOR LAYER-MEAN TRANSMISSIVITIES.
!
!          FOR A GIVEN CO2 MIXING RATIO AND STANDARD TEMPERATURE
!      PROFILE,A TABLE OF TRANSMISSION FUNCTIONS FOR A FIXED GRID
!     OF ATMOSPHERIC PRESSURES HAS BEEN PRE-CALCULATED.
!      THE STANDARD TEMPERATURE PROFILE IS COMPUTED FROM THE US
!     STANDARD ATMOSPHERE (1977) TABLE.ADDITIONALLY, THE
!     SAME TRANSMISSION FUNCTIONS HAVE BEEN PRE-CALCULATED FOR A
!     TEMPERATURE PROFILE INCREASED AND DECREASED (AT ALL LEVELS)
!     BY 25 DEGREES.
!         THIS PROGRAM READS IN THE PRESPECIFIED TRANSMISSION FUNCTIONS
!     AND A USER-SUPPLIED PRESSURE GRID (P(I)) AND CALCULATES TRANS-
!     MISSION FUNCTIONS ,TAU(P(I),P(J)), FOR ALL P(I) S AND P(J) S.
!     A LOGARITHMIC INTERPOLATION SCHEME IS USED.
!         THIS METHOD IS REPEATED FOR THE THREE TEMPERATURE PROFILES
!     GIVEN ABOVE .THEREFORE OUTPUTS FROM THE PROGRAM ARE THREE TABLES
!     OF TRANSMISSION FUNCTIONS FOR THE USER-SUPPLIED PRESSURE GRID.
!     THE EXISTENCE OF THE THREE TABLES PERMITS SUBSEQUENT INTERPO-
!     LATION TO A USER-SUPPLIED TEMPERATURE PROFILE USING THE METHOD
!     DESCRIBED IN THE REFERENCE.SEE LIMITATIONS SECTION IF THE
!     USER DESIRES TO OBTAIN ONLY 1 TABLE OF TRANSMISSIVITIES.
!
!     MODIFICATIONS FOR LAYER-MEAN TRANSMISSIVITIES:
!          THE PRESSURES INPUTTED ARE THE LAYER-MEAN PRESSURES,PD,
!     AND THE LAYER-MEAN PRESSURES ,PLM. A SERIES OF TRANSMISSIVITIES
!     (TAU(P'',PLM(J)) ARE COMPUTED AND THE INTEGRAL GIVEN IN THE
!     DISCUSSION OF PURPOSE 2 IS COMPUTED.FOR PLM(I) NOT EQUAL TO
!     PLM(J) SIMPSON S RULE IS USED WITH 5 POINTS. IF PLM(I)=PLM(J)
!     (THE -NEARBY LAYER- CASE) A 49-POINT QUADRATURE IS USED FOR
!     GREATER ACCURACY.THE OUTPUT IS IN TAULM(PLM(I),PLM(J)).
!        NOTE:
!     TAULM IS NOT A SYMMETRICAL MATRIX. FOR THE ARRAY ELEMENT
!     TAULM(PLM(I),PLM(J)),THE INNER(FIRST,MOST RAPIDLY VARYING)
!     DIMENSION IS THE VARYING LAYER-MEAN PRESSURE,PLM(I);THE OUTER
!     (SECOND) DIMENSION IS THE FIXED LAYER-MEAN PRESSURE PLM(J).
!     THUS THE ELEMENT TAULM(2,3) IS THE TRANSMISSION FUNCTION BETWEEN
!     THE FIXED PRESSURE PLM(3)  AND THE PRESSURE LAYER HAVING AN AVERAG
!     PRESSURE OF PLM(2).
!         ALSO NOTE THAT NO QUADRATURE IS PERFORMED OVER THE LAYER
!     BETWEEN THE SMALLEST NONZERO PRESSURE AND ZERO PRESSURE;
!     TAULM IS TAULM(0,PLM(J)) IN THIS CASE,AND TAULM(0,0)=1.
!
!
!             REFERENCE:
!         S.B.FELS AND M.D.SCHWARZKOPF,-AN EFFICIENT ACCURATE
!     ALGORITHM FOR CALCULATING CO2 15 UM BAND COOLING RATES-,JOURNAL
!     OF GEOPHYSICAL RESEARCH,VOL.86,NO. C2, PP.1205-1232,1981.
!        MODIFICATIONS TO THE ALGORITHM HAVE BEEN MADE BY THE AUTHORS;
!     CONTACT S.B.F.OR M.D.S. FOR FURTHER DETAILS.A NOTE TO J.G.R.
!     IS PLANNED TO DOCUMENT THESE CHANGES.
!
!            AUTHOR:    M.DANIEL SCHWARZKOPF
!
!            DATE:      14 JULY 1983
!
!            ADDRESS:
!
!                      G.F.D.L.
!                      P.O.BOX 308
!                      PRINCETON,N.J.08540
!                      U.S.A.
!            TELEPHONE:  (609) 452-6521
!
!            INFORMATION ON TAPE: THIS SOURCE IS THE FIRST FILE
!        ON THIS TAPE.THE SIX FILES THAT FOLLOW ARE CO2 TRANS-
!        MISSIVITIES FOR THE 500-850 CM-1 INTERVAL FOR CO2
!        CONCENTRATIONS OF 330 PPMV (1X) ,660 PPMV (2X), AND
!        1320 PPMV (4X). THE FILES ARE ARRANGED AS FOLLOWS:
!          FILE 2   1X,CONSOLIDATED USING B(250) WEIGHTING FCTN.
!          FILE 3   1X,CONSOLIDATED WITH NO WEIGHTING FCTN.
!          FILE 4   2X,CONSOLIDATED USING B(250) WEIGHTING FCTN.
!          FILE 5   2X,CONSOLIDATED WITH NO WEIGHTING FCTN.
!          FILE 6   4X,CONSOLIDATED USING B(250) WEIGHTING FCTN.
!          FILE 7   4X,CONSOLIDATED WITH NO WEIGHTING FCTN.
!            FILES 2,4,6 ARE RECOMMENDED FOR USE IN OBTAINING
!        TRANSMISSION FUNCTIONS FOR USE IN HEATING RATE
!        COMPUTATIONS;THEY CORRESPOND TO THE TRANSMISSIVITIES
!        DISCUSSED IN THE 1980 PAPER.FILES 3,5,7 ARE PROVIDED
!        TO FACILITATE COMPARISON WITH OBSERVATION AND WITH OTHER
!        CALCULATIONS.
!
!            PROGRAM LANGUAGE: FORTRAN 1977,INCLUDING PARAMETER
!        AND PROGRAM STATEMENTS.THE PROGRAM IS WRITTEN ON A
!        CYBER 170-730.SEE THE SECTION ON LIMITATIONS FOR
!        ADAPTATIONS TO OTHER MACHINES.
!
!          INPUT UNITS,FORMATS AND FORMAT STATEMENT NOS:
!
!   UNIT NO    VARIABLES       FORMAT      STATEMENT NO.    TYPE
!      5        P (PURPOSE 1)  (5E16.9)        201         CARDS
!      5        PD (PURPOSE 2) (5E16.9)        201         CARDS
!      5        PLM(PURPOSE 2) (5E16.9)        201         CARDS
!      5        NMETHD         (I3)            202         CARDS
!      20       TRANSA         (4F20.14)       102          TAPE
!NOV89
!      ITAPE    TRANSA         (4F20.14)       102          TAPE
!NOV89
!
!         OUTPUT UNITS,FORMATS AND FORMAT STATEMENT NOS:
!
!   UNIT NO    VARIABLES       FORMAT     STATEMENT NO.
!      6         TRNFCT        (1X,8F15.8)     301         PRINT
!      22        TRNFCT        (4F20.14)       102          TAPE
!
!            PARAMETER INPUTS:
!     A) NLEVLS    : NLEVLS IS AN (INTEGER) PARAMETER DENOTING
!        THE NUMBER OF NONZERO PRESSURE LEVELS FOR PURPOSE 1
!        OR THE NUMBER OF NONZERO LAYER PRESSURES NEEDED TO
!        SPECIFY THE PRESSURE LAYERS(PURPOSE 2) IN THE OUTPUT
!        GRID. FOR EXAMPLE,IN PURPOSE 1,IF P=0,100,1000,NLEVLS=2.
!        IF,IN PURPOSE 2,PD=0,100,500,1000,THE NUMBER OF NONZERO
!        PRESSURE LAYERS=2,SO NLEVLS=2
!           IN THE CODE AS WRITTEN,NLEVLS=40; THE USER SHOULD
!        CHANGE THIS VALUE TO A USER-SPECIFIED VALUE.
!     B) NLP1,NLP2 : INTEGER PARAMETERS DEFINED AS: NLP1=NLEVLS+1;
!        NLP2=NLEVLS+2.
!           SEE LIMITATIONS FOR CODE MODIFICATIONS IF PARAMETER
!        STATEMENTS ARE NOT ALLOWED ON YOUR MACHINE.
!
!            INPUTS:
!
!     A) TRANSA    : THE 109X109 GRID OF TRANSMISSION FUNCTIONS
!            TRANSA IS A  DOUBLE PRECISION REAL ARRAY.
!
!           TRANSA  IS READ FROM FILE 20. THIS FILE CONTAINS 3
!     RECORDS,AS FOLLOWS:
!        1)   TRANSA, STANDARD TEMPERATURE PROFILE
!        3)   TRANSA, STANDARD TEMPERATURES + 25 DEG
!        5)   TRANSA, STANDARD TEMPERATURES - 25 DEG
!
!     B)   NMETHD: AN INTEGER WHOSE VALUE IS EITHER 1 (IF CO2INT IS
!       TO BE USED FOR PURPOSE 1) OR 2 (IF CO2INT IS TO BE USED FOR
!       PURPOSE 2).
!
!     C)     P,PD,PLM :
!          P IS A REAL ARRAY (LENGTH NLP1) SPECIFYING THE PRESSURE
!       GRID AT WHICH TRANSMISSION FUNCTIONS ARE TO BE COMPUTED FOR
!       PURPOSE 1.THE DIMENSION  OF P IS  IN MILLIBARS.THE
!       FOLLOWING LIMITATIONS WILL BE EXPLAINED MORE
!       IN THE SECTION ON LIMITATIONS: P(1) MUST BE ZERO; P(NLP1),THE
!       LARGEST PRESSURE, MUST NOT EXCEED 1165 MILLIBARS.
!         PD IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE PRESSURE
!       LAYERS FOR WHICH LAYER-AVERAGED TRANSMISSION FUNCTIONS ARE
!       TO BE COMPUTED.THE DIMENSION OF PD IS MILLIBARS.THE LIMITATIONS
!       FOR PD ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON
!       LIMITATIONS.
!         PLM IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE LAYER-MEAN
!       PRESSURES. THE DIMENSION OF PLM IS MILLIBARS. THE LIMITATIONS
!       FOR PLM ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON
!       LIMITATIONS.PD IS READ IN BEFORE PLM.
!
!          NOTE: AGAIN,WE NOTE THAT THE USER WILL INPUT EITHER P (FOR
!       PURPOSE 1) OR PD AND PLM(FOR PURPOSE 2) BUT NOT BOTH.
!
!
!
!
!           LIMITATIONS:
!     1)       P(1)=0.,PD(1)=0.,PLM(1)=0. THE TOP PRESSURE LEVEL
!       MUST BE ZERO,OR THE TOP PRESSURE LAYER MUST BE BOUNDED BY ZERO.
!       THE TOP LAYER-MEAN PRESSURE (PLM(1)) MUST BE ZERO; NO
!       QUADRATURE IS DONE ON THE TOP PRESSURE LAYER.EVEN IF ONE IS
!       NOT INTERESTED IN THE TRANSMISSION FUNCTION BETWEEN 0 AND P(J),
!       ONE MUST INCLUDE SUCH A LEVEL.
!     2)      PD(NLP2)=P(NLP1) IS LESS THAN OR EQUAL TO 1165 MB.
!       EXTRAPOLATION TO HIGHER PRESSURES IS NOT POSSIBLE.
!     3)      IF PROGRAM IS NOT PERMITTED ON YOUR COMPILER,
!       SIMPLY DELETE THE LINE.
!     4)      IF PARAMETER IS NOT PERMITTED,DO THE FOLLOWING:
!            1) DELETE ALL PARAMETER STATEMENTS IN CO2INT
!            2) AT THE POINT WHERE NMETHOD IS READ IN,ADD:
!                READ (5,202) NLEVLS
!                NLP1=NLEVLS+1
!                NLP2=NLEVLS+2
!            3) CHANGE DIMENSION AND/OR COMMON STATEMENTS DEFINING
!              ARRAYS TRNS,DELTA,P,PD,TRNFCT,PS,PDS,PLM IN CO2INT.
!              THE NUMERICAL VALUE OF (NLEVLS+1) SHOULD BE INSERTED
!              IN DIMENSION OR COMMON STATEMENTS FOR TRNS,DELTA,
!              P,TRNFCT,PS,PLM; THE NUMERICAL VALUE OF (NLEVLS+2)
!              IN DIMENSION OR COMMON STATEMENTS FOR PD,PDS.
!      5)    PARAMETER (NLEVLS=40) AND THE OTHER PARAMETER
!       STATEMENTS ARE WRITTEN IN CDC FORTRAN; ON OTHER MACHINES THE
!       SAME STATEMENT MAY BE WRITTEN DIFFERENTLY,FOR EXAMPLE AS
!       PARAMETER   NLEVLS=40
!      6) -DOUBLE PRECISION- IS USED INSTEAD OF -REAL*8- ,DUE TO
!       REQUIREMENTS OF CDC FORTAN.
!      7) THE STATEMENT -DO 400 KKK=1,3- CONTROLS THE NUMBER OF
!       TRANSMISSIVITY OUTPUT MATRICES PORDUCED BY THE PROGRAM.TO
!       PRODUCE 1 OUTPUT MATRIX,DELETE THIS STATEMENT.
!
!     OUTPUT:
!         A) TRNFCT IS AN (NLP1,NLP1) REAL ARRAY OF THE TRANSMISSION
!     FUNCTIONS APPROPRIATE TO YOUR ARRAY. IT IS TO BE SAVED ON FILE 22.
!     THE PROCEDURE FOR SAVING MAY BE MODIFIED; AS GIVEN HERE,THE
!     OUTPUT IS IN CARD IMAGE FORM WITH A FORMAT OF (4F20.14).
!
!         B)  PRINTED  OUTPUT IS A LISTING OF TRNFCT ON UNIT 6, IN
!     THE FORMAT (1X,8F15.8) (FORMAT STATEMENT 301). THE USER MAY
!     MODIFY OR ELIMINATE THIS AT WILL.
!
!      ************   FUNCTION INTERPOLATER ROUTINE  *****************
!
!
!     ******   THE FOLLOWING PARAMETER GIVES THE NUMBER OF     *******
!     ******           DATA LEVELS IN THE MODEL                *******
!     ****************************************************************
!     ****************************************************************
      COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N
!     COMMON/PRESS/PA(109)
!     COMMON/TRAN/ TRANSA(109,109)
!     COMMON / OUTPUT / TRNS(NLP1,NLP1)
!     COMMON/INPUTP/P(NLP1),PD(NLP2)
      DIMENSION TRNS(NLP1,NLP1)
      DIMENSION P(NLP1),PD(NLP2)
      DIMENSION PS(NLP1),PDS(NLP2),PLM(NLP1)
      DIMENSION NRTAB(3)
      DIMENSION T15A(NLP2,2),T15B(NLP1)
      DIMENSION T22(NLP1,NLP1,3)
      LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
      DATA NRTAB/1,2,4/
!***********************************
!   THE FOLLOWING ARE THE INPUT FORMATS
100   FORMAT (4F20.14)
743   FORMAT (F20.14)
201   FORMAT (5E16.9)
202   FORMAT (I3)
!O222   203   FORMAT (F12.6,I2)
203   FORMAT (F12.6)
!    THE FOLLOWING ARE THE OUTPUT FORMATS
102   FORMAT (4F20.14)
301   FORMAT (1X,8F15.8)
!
!CC   REWIND 15
!CC   REWIND 20
!NOV89
      REWIND ITAPE
!NOV89
!CC   REWIND 22
!
!     CALCULATION OF PA -THE -TABLE- OF 109 GRID PRESSURES
!     NOTE-THIS CODE MUST NOT BE CHANGED BY THE USER^^^^^^^^^
      PA(1)=0.
      FACT15=10.**(1./15.)
      FACT30=10.**(1./30.)
      PA(2)=1.0E-3
      DO 231 I=2,76
      PA(I+1)=PA(I)*FACT15
231   CONTINUE
      DO 232 I=77,108
      PA(I+1)=PA(I)*FACT30
232   CONTINUE
!
      N=25
      NLV=NLEVLS
      NLP1V=NLP1
      NLP2V=NLP2
!     READ IN THE CO2 MIXING RATIO(IN UNITS OF 330 PPMV),AND AN INDEX
!     GIVING THE FREQUENCY RANGE OF THE LBL DATA
!O222    READ (5,203) RATIO,IR
!CC         IR = 1
!CC         READ (5,203) RATIO
!O222   ***********************************
!***VALUES FOR IR*****
!          IR=1     CONSOL. LBL TRANS. =490-850
!          IR=2     CONSOL. LBL TRANS. =490-670
!          IR=3     CONSOL. LBL TRANS. =670-850
!          IR=4     CONSOL. LBL TRANS. =2270-2380
!*** IR MUST BE 1,2,3 OR 4 FOR THE PGM. TO WORK
!     ALSO READ IN THE METHOD NO.(1 OR 2)
!CC         READ (5,202) NMETHD
      GO TO 621
!NOV89  621   ITAP1=20
621   ITAP1=ITAPE
!NOV89
      NTAP=1
      IF (NMETHD.EQ.2) GO TO 502
!   *****CARDS FOR PURPOSE 1(NMETHD=1)
!CC         READ (15,201) (P(I),I=1,NLP1)
      DO 300 I=1,NLP1
        P(I)=T15B(I)
  300 CONTINUE
      DO 801 I=1,NLP1
      PS(I)=P(I)
801   CONTINUE
      GO TO 503
502   CONTINUE
!  *****CARDS FOR PURPOSE 2(NMETHD=2)
!CC         READ (15,201) (PD(I),I=1,NLP2)
!CC         READ (15,201) (PLM(I),I=1,NLP1)
      DO 303 I=1,NLP2
        PD(I)=T15A(I,1)
  303 CONTINUE
      DO 302 I=1,NLP1
        PLM(I)=T15A(I,2)
  302 CONTINUE
      DO 802 I=1,NLP1
      PDS(I)=PD(I+1)
      PS(I)=PLM(I)
802   CONTINUE
!
503   CONTINUE
!  *****DO LOOP CONTROLLING NUMBER OF OUTPUT MATRICES
!NOV89
!NOV89    DO 400 KKK=1,3
      ICLOOP = 3
      IF (IR.EQ.4) ICLOOP = 1
      DO 400 KKK=1,ICLOOP
!NOV89
!  **********************
      IF (NMETHD.EQ.2) GO TO 505
!   *****CARDS FOR PURPOSE 1(NMETHD=1)
      DO 803 I=1,NLP1
      P(I)=PS(I)
803   CONTINUE
      GO TO 506
505   CONTINUE
!  *****CARDS FOR PURPOSE 2(NMETHD=2)
      DO 804 I=1,NLP1
      PD(I)=PDS(I)
      P(I)=PS(I)
804   CONTINUE
!
506   CONTINUE
      IA=108
      IAP=IA+1
!NOV89   IF (NTAP.EQ.1) READ (20,100) ((TRANSA(I,J),I=1,109),J=1,109)
!mp       IF (NTAP.EQ.1) READ (ITAPE,100) ((TRANSA(I,J),I=1,109),J=1,109)
        IF (NTAP.EQ.1) THEN
           IF ( wrf_dm_on_monitor() ) READ (ITAPE,743) ((TRANSA(I,J),I=1,109),J=1,109)
           CALL wrf_dm_bcast_bytes ( TRANSA , size ( TRANSA ) * RWORDSIZE )
        ENDIF
!mp	IF (NTAP.EQ.1) READ (ITAPE,100) (tmp(I),I=1,11881
!mp
	do J=109,1,-6
!mp	write(6,697)(TRANSA(I,J),I=5,105,10)
	enddo
! 697	format(11(f5.3,1x))
!mp
!NOV89
      DO 4 I=1,IAP
      TRANSA(I,I)=1.0
    4 CONTINUE
      CALL COEINT(RATIO,IR)
      DO 805 I=1,NLP1
      DO 805 J=1,NLP1
      TRNS(J,I)=1.00
805   CONTINUE
      DO 10 I=1,NLP1
      DO 20 J=1,I
      IF (I.EQ.J) GO TO 20
      P1=P(J)
      P2=P(I)
      CALL SINTR2
      TRNS(J,I)=TRNSLO
20    CONTINUE
10    CONTINUE
      DO 47 I=1,NLP1
      DO 47 J=I,NLP1
      TRNS(J,I)=TRNS(I,J)
47    CONTINUE
!  *****THIS IS THE END OF PURPOSE 1 CALCULATIONS
      IF (NMETHD.EQ.1) GO TO 2872
!
      DO 51 J=1,NLP1
      DO 52 I=2,NLP1
      IA=I
      JA=J
      N=25
      IF (I.NE.J) N=3
      CALL QUADSR(NLV,NLP1V,NLP2V,P,PD,TRNS)
52    CONTINUE
51    CONTINUE
!  *****THIS IS THE END OF PURPOSE 2 CALCULATIONS
2872  CONTINUE
!
!+++  WRITE (6,301) ((TRNS(I,J),I=1,NLP1),J=1,NLP1)
!CC         WRITE (22,102) ((TRNS(I,J),I=1,NLP1),J=1,NLP1)
      DO 304 J=1,NLP1
       DO 304 I=1,NLP1
        T22(I,J,KKK) = TRNS(I,J)
  304 CONTINUE
400   CONTINUE
      RETURN
      END SUBROUTINE CO2INT
!CCC  PROGRAM CO2IN1
      SUBROUTINE CO2IN1(T20,T21,T66,IQ,L,LP1)
!    CO2IN1=CO2INS FOR METHOD 1
!     *********************************************************
!       SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 ***
!          ..... K.CAMPANA   MARCH 1988,OCTOBER 1988
!          ..... K.CAMPANA   DECEMBER 88 CLEANED UP FOR LAUNCHER
!     *********************************************************
      DIMENSION T20(LP1,LP1,3),T21(LP1,LP1,3),T66(L,6)
      DIMENSION DCDT8(LP1,LP1),DCDT10(LP1,LP1),CO2PO(LP1,LP1), &
       CO2800(LP1,LP1),CO2PO1(LP1,LP1),CO2801(LP1,LP1),CO2PO2(LP1,LP1), &
       CO2802(LP1,LP1),N(LP1),D2CT8(LP1,LP1),D2CT10(LP1,LP1)
      ITIN=20
      ITIN1=21
!O222 LATEST CODE HAS IQ=1
!CC         IQ=4
1011  FORMAT (4F20.14)
!CC        READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1)
!CC        READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1)
!CC        READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1)
!CC        READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1)
!CC        READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1)
!CC        READ (ITIN1,1011) ((CO2802(I,J),I=1,LP1),J=1,LP1)
      DO 300 J=1,LP1
        DO 300 I=1,LP1
          CO2PO(I,J) = T20(I,J,1)
!NOV89
          IF (IQ.EQ.5) GO TO 300
!NOV89
          CO2PO1(I,J) = T20(I,J,2)
          CO2PO2(I,J) = T20(I,J,3)
  300 CONTINUE
      DO 301 J=1,LP1
        DO 301 I=1,LP1
          CO2800(I,J) = T21(I,J,1)
!NOV89
          IF (IQ.EQ.5) GO TO 301
!NOV89
          CO2801(I,J) = T21(I,J,2)
          CO2802(I,J) = T21(I,J,3)
  301 CONTINUE
!***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS
!   ARE:
!        IQ=1    560-800     (CONSOL.=490-850)
!        IQ=2    560-670     (CONSOL.=490-670)
!        IQ=3    670-800     (CONSOL.=670-850)
!        IQ=4    560-760 (ORIGINAL CODE)   (CONSOL.=490-850)
!NOV89
!        IQ=5   2270-2380    (CONSOL.=2270-2380)
!NOV89
!  THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS
!  USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT
!  WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S.
      IF (IQ.EQ.1) THEN
         C1=1.5
         C2x=0.5
      ENDIF
      IF (IQ.EQ.2) THEN
        C1=18./11.
        C2x=7./11.
      ENDIF
      IF (IQ.EQ.3) THEN
        C1=18./13.
        C2x=5./13.
      ENDIF
      IF (IQ.EQ.4) THEN
        C1=1.8
        C2x=0.8
      ENDIF
!NOV89
      IF (IQ.EQ.5) THEN
        C1=1.0
        C2x=0.0
      ENDIF
!NOV89
      DO 1021 I=1,LP1
      DO 1021 J=1,LP1
      CO2PO(J,I)=C1*CO2PO(J,I)-C2x
      CO2800(J,I)=C1*CO2800(J,I)-C2x
!NOV89
      IF (IQ.EQ.5) GO TO 1021
!NOV89
      CO2PO1(J,I)=C1*CO2PO1(J,I)-C2x
      CO2801(J,I)=C1*CO2801(J,I)-C2x
      CO2PO2(J,I)=C1*CO2PO2(J,I)-C2x
      CO2802(J,I)=C1*CO2802(J,I)-C2x
1021  CONTINUE
!NOV89
      IF (IQ.GE.1.AND.IQ.LE.4) THEN
!NOV89
      DO 1 J=1,LP1
      DO 1 I=1,LP1
      DCDT8(I,J)=.02*(CO2801(I,J)-CO2802(I,J))*100.
      DCDT10(I,J)=.02*(CO2PO1(I,J)-CO2PO2(I,J))*100.
      D2CT8(I,J)=.0016*(CO2801(I,J)+CO2802(I,J)-2.*CO2800(I,J))*1000.
      D2CT10(I,J)=.0016*(CO2PO1(I,J)+CO2PO2(I,J)-2.*CO2PO(I,J))*1000.
1     CONTINUE
!NOV89
      ENDIF
!NOV89
!O222 *********************************************************
!CC          REWIND 66
!        SAVE CDTM51,CO2M51,C2DM51,CDTM58,CO2M58,C2DM58..ON TEMPO FILE
!CC          WRITE (66) (DCDT10(I,I+1),I=1,L)
!CC          WRITE (66) (CO2PO(I,I+1),I=1,L)
!CC          WRITE (66) (D2CT10(I,I+1),I=1,L)
!CC          WRITE (66) (DCDT8(I,I+1),I=1,L)
!CC          WRITE (66) (CO2800(I,I+1),I=1,L)
!CC          WRITE (66) (D2CT8(I,I+1),I=1,L)
!CC          REWIND 66
!O222 *********************************************************
      DO 400 I=1,L
        T66(I,2) = CO2PO(I,I+1)
        T66(I,5) = CO2800(I,I+1)
!NOV89
        IF (IQ.EQ.5) GO TO 400
!NOV89
        T66(I,1) = DCDT10(I,I+1)
        T66(I,3) = D2CT10(I,I+1)
        T66(I,4) = DCDT8(I,I+1)
        T66(I,6) = D2CT8(I,I+1)
  400 CONTINUE
      RETURN
      END SUBROUTINE CO2IN1
!CCC  PROGRAM PTZ - COURTESY OF DAN SCHWARZKOPF,GFDL DEC 1987....
      SUBROUTINE CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, &
                        SFULL,SHALF,PPTOP,LREAD,NL,NLP,NLP2)
!
! **         THIS PROGRAM CALCULATES TEMPERATURES ,H2O MIXING RATIOS
! **         AND O3 MIXING RATIOS BY USING AN ANALYTICAL
! **         FUNCTION WHICH APPROXIMATES
! **         THE US STANDARD (1976).  THIS IS
! **         CALCULATED IN FUNCTION 'ANTEMP', WHICH IS CALLED BY THE
! **         MAIN PROGRAM.  THE FORM OF THE ANALYTICAL FUNCTION WAS
! **         SUGGESTED TO ME IN 1971 BY RICHARD S. LINDZEN.
! ******************************************************************
!         CODE TO SAVE STEMP,GTEMP ON DATA SET,BRACKETED BY CO222  **
!             ....K. CAMPANA MARCH 88,OCTOBER 88
      DIMENSION SGTEMP(NLP,2),T41(NLP2,2),T42(NLP), &
                T43(NLP2,2),T44(NLP)
      DIMENSION SGLVNU(NLP),SIGLNU(NL)
      DIMENSION SFULL(NLP),SHALF(NL)
! ******************************************************************
!
!*****THIS VERSION IS ONLY USABLE FOR 1976 US STD ATM AND OBTAINS
!     QUANTITIES FOR CO2 INTERPOLATION AND INSERTION INTO OPERA-
!     TIONAL RADIATION CODES
!
      CHARACTER*20 PROFIL
      DIMENSION PRESS(NLP),TEMP(NLP),ALT(NLP),WMIX(NLP),O3MIX(NLP)
      DIMENSION WMXINT(NLP,4),WMXOUT(NLP2),OMXINT(NLP,4),OMXOUT(NLP2)
      DIMENSION PD(NLP2),GTEMP(NLP)
      DIMENSION PRS(NLP),TEMPS(NLP),PRSINT(NLP),TMPINT(NLP,4),A(NLP,4)
      DIMENSION PROUT(NLP2),TMPOUT(NLP2),TMPFLX(NLP2),TMPMID(NLP2)
!
!
      DATA PROFIL/ &
         'US STANDARD 1976'/
      DATA PSMAX/1013.250/
!
! **         NTYPE IS AN INTEGER VARIABLE WHICH HAS THE FOLLOWING
! **        VALUES:    0 =SIGMA LEVELS ARE USED;   1= SKYHI L40 LEVELS
! **        ARE USED;   2 = SKYHI L80 LEVELS ARE USED. DEFAULT: 0
!
      NTYPE=0
!O222 READ (*,*) NTYPE
    5 NLEV=NL
      DELZAP=0.5
      R=8.31432
      G0=9.80665
      ZMASS=28.9644
      AA=6356.766
         ALT(1)=0.0
         TEMP(1)=ANTEMP(6,0.0)
!*******DETERMINE THE PRESSURES (PRESS)
      PSTAR=PSMAX
!
!***  LTOP COMPUTATION MOVED FROM MODEL INITIALIZATION
!
      LTOP(1)=0
      LTOP(2)=0
      LTOP(3)=0
      DO 30 N=1,NL
        PCLD=(PSTAR-PPTOP*10.)*SHALF(N)+PPTOP*10.
        IF(PCLD.GE.642.)LTOP(1)=N
        IF(PCLD.GE.350.)LTOP(2)=N
        IF(PCLD.GE.150.)LTOP(3)=N
!       PRINT *,N,PCLD,SHALF(N),PSTAR,PPTOP
   30 CONTINUE
!
!O222 IF (NTYPE.EQ.1) CALL SKYP(PSTAR,PD,GTEMP)
!O222 IF (NTYPE.EQ.2) CALL SKY80P(PSTAR,PD,GTEMP)
!O222 IF (NTYPE.EQ.0) CALL SIGP(PSTAR,PD,GTEMP)
!CC----      CALL SIGP(PSTAR,PD,GTEMP)
      NLM=NL-1
      CALL SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, &
                SFULL,SHALF,PPTOP,LREAD,NL,NLP,NLM,NLP2)
      PD(NLP2)=PSTAR
      DO 40 N=1,NLP
      PRSINT(N)=PD(NLP2+1-N)
 40   CONTINUE
!    *** CALCULATE TEMPS FOR SEVERAL PRESSURES TO DO QUADRATURE
      DO 504 NQ=1,4
      DO 505 N=2,NLP
 505  PRESS(N)=PRSINT(N)+0.25*(NQ-1)*(PRSINT(N-1)-PRSINT(N))
      PRESS(1)=PRSINT(1)
!*********************
      DO 100 N=1,NLEV
!
! **         ESTABLISH COMPUTATATIONAL LEVELS BETWEEN USER LEVELS AT
! **         INTERVALS OF APPROXIMATELY 'DELZAP' KM.
!
      DLOGP=7.0*ALOG(PRESS(N)/PRESS(N+1))
      NINT=DLOGP/DELZAP
      NINT=NINT+1
      ZNINT=NINT
!     G=G0
      DZ=R*DLOGP/(7.0*ZMASS*G0*ZNINT)
      HT=ALT(N)
!
! **         CALCULATE HEIGHT AT NEXT USER LEVEL BY MEANS OF
! **                   RUNGE-KUTTA INTEGRATION.
!
      DO 200 M=1,NINT
      RK1=ANTEMP(6,HT)*DZ
      RK2=ANTEMP(6,HT+0.5*RK1)*DZ
      RK3=ANTEMP(6,HT+0.5*RK2)*DZ
      RK4=ANTEMP(6,HT+RK3)*DZ
!mp	write(6,*) 'RK values,DZ ', RK1,RK2,RK3,RK4,DZ
      HT=HT+0.16666667*(RK1+RK2+RK2+RK3+RK3+RK4)
  200 CONTINUE
      ALT(N+1)=HT
      TEMP(N+1)=ANTEMP(6,HT)
  100 CONTINUE
      DO 506 N=1,NLP
      TMPINT(N,NQ)=TEMP(N)
      A(N,NQ)=ALT(N)
506   CONTINUE
504   CONTINUE
!O222   *****************************************************
!***OUTPUT TEMPERATURES
!O222   *****************************************************
      DO 901 N=1,NLP
        SGTEMP(N,1) = TMPINT(NLP2-N,1)
  901 CONTINUE
!O222   *****************************************************
!***OUTPUT GTEMP
!O222   *****************************************************
      DO 902 N=1,NLP
        SGTEMP(N,2) = GTEMP(N)
  902 CONTINUE
!O222   *****************************************************
      RETURN
      END SUBROUTINE CO2PTZ
      FUNCTION PATH(A,B,C,E)
!....
!     DOUBLE PRECISION XA,CA
!     COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP
      PEXP=1./SEXP
      PATH=((A-B)**PEXP*(A+B+C))/(E*(A+B+C)+(A-B)**(PEXP-1.))
      RETURN
      END FUNCTION PATH
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE QINTRP(XM,X0,XP,FM,F0,FP,X,F)
!....
!     DOUBLE PRECISION FM,F0,FP,F,D1,D2,B,A,DEL
      D1=(FP-F0)/(XP-X0)
      D2=(FM-F0)/(XM-X0)
      B=(D1-D2)/(XP-XM)
      A=D1-B*(XP-X0)
      DEL=(X-X0)
      F=F0+DEL*(A+DEL*B)
      RETURN
      END SUBROUTINE QINTRP
      SUBROUTINE QUADSR(NLV,NLP1V,NLP2V,P,PD,TRNS)
      COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N
      DIMENSION P(NLP1V),PD(NLP2V),TRNS(NLP1V,NLP1V)
      DIMENSION WT(101)
      N2=2*N
      N2P=2*N+1
!  *****WEIGHTS ARE CALCULATED
      WT(1)=1.
      DO 21 I=1,N
      WT(2*I)=4.
      WT(2*I+1)=1.
21    CONTINUE
      IF (N.EQ.1) GO TO 25
      DO 22 I=2,N
      WT(2*I-1)=2.
22    CONTINUE
25    CONTINUE
      TRNSNB=0.
      DP=(PD(IA)-PD(IA-1))/N2
      PFIX=P(JA)
      DO 1 KK=1,N2P
      PVARY=PD(IA-1)+(KK-1)*DP
      IF (PVARY.GE.PFIX) P2=PVARY
      IF (PVARY.GE.PFIX) P1=PFIX
      IF (PVARY.LT.PFIX) P1=PVARY
      IF (PVARY.LT.PFIX) P2=PFIX
      CALL SINTR2
      TRNSNB=TRNSNB+TRNSLO*WT(KK)
1     CONTINUE
      TRNS(IA,JA)=TRNSNB*DP/(3.*(PD(IA)-PD(IA-1)))
      RETURN
      END SUBROUTINE QUADSR
!---------------------------------------------------------------------
      SUBROUTINE SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, &
                      SIGLV,SIGLY,PPTOP,LREAD,KD,KP,KM,KP2)
      DIMENSION Q(KD),QMH(KP),PD(KP2),PLM(KP),GTEMP(KP),PDT(KP2)
      DIMENSION SIGLY(KD),SIGLV(KP)
      DIMENSION CI(KP),SGLVNU(KP),DEL(KD),SIGLNU(KD),CL(KD),RPI(KM)
      DIMENSION IDATE(4)
      DIMENSION T41(KP2,2),T42(KP), &
                T43(KP2,2),T44(KP)
!     integer :: retval
!     character(50) :: prsmid='prsmid'
!CC   18 LEVEL SIGMAS FOR NMC MRF(NEW) MODEL
!CC   DATA Q/.021,.074,.124,.175,.225,.275,.325,.375,.425,.497, &
!CC          .594,.688,.777,.856,.920,.960,.981,.995/
!     FOR SIGMA MODELS,Q=SIGMA,QMH=0.5(Q(I)+Q(I+1),
!     PD=Q*PSS,PLM=QMH*PSS.PSS=SURFACE PRESSURE(SPEC.)
!
!.....   GET NMC SIGMA STRUCTURE
!CC   IF (LREAD.GT.0) GO TO 914
!---   PPTOP IS MODEL TOP PRESSURE IN CB....
!        SIGMA DATA IS BOTTOM OF ATMOSPHERE TO T.O.A.....
!cccc PPTOP=5.0
!     READ(11,PPTOP,END=12321)
12321 CONTINUE
!     WRITE(6,88221)PPTOP,KD,KP
!88221 FORMAT(' ENTER SIGP PPTOP=',E12.5,' KD=',I2,' KP=',I2)
!     open(unit=23,file='fort.23',form='unformatted' &
!     ,    access='sequential')
!     REWIND 23
!     READ(23)SIGLY
!     DO KKK=1,KD
!      SIGLY(KKK)=1.-(FLOAT(KKK)-0.5)/KD
!     END DO
!     WRITE(6,88222)
!88222 FORMAT(' READ AETA')
!     DO 37821 LLL=1,KD
!     WRITE(6,37820)LLL,SIGLY(LLL)
!37820 FORMAT(' L=',I2,' AETA=',E12.5)
!37821 CONTINUE
!     READ(23)SIGLV
!     DO KKK=1,KP
!      SIGLV(KKK)=1.-(FLOAT(KKK-1))/KD
!     END DO
!     WRITE(6,88223)
!88223 FORMAT(' READ ETA')
!     PRINT 704,(SIGLY(K),K=1,KD)
!     PRINT 704,(SIGLV(K),K=1,KP)
!      DO 37823 LLL=1,KP
!      WRITE(6,37822)LLL,SIGLV(LLL)
!37822 FORMAT(' L=',I2,' ETA=',E12.5)
!37823 CONTINUE
  701 FORMAT(F6.2)
  702 FORMAT(7F10.6)
      IF (PPTOP.LE.0.) GO TO 708
      PSFC=100.
!--- IF PTOP NOT EQUAL TO ZERO ADJUST SIGMA SO AS TO GET PROPER STD ATM
!      VERTICAL LOCATION
      DO 706 K=1,KD
       SIGLY(K) = (SIGLY(K)*(PSFC-PPTOP)+PPTOP)/PSFC
  706 CONTINUE
      DO 707 K=1,KP
       SIGLV(K) = (SIGLV(K)*(PSFC-PPTOP)+PPTOP)/PSFC
  707 CONTINUE
  708 CONTINUE
!     PRINT 703,PPTOP
!     PRINT 704,(SIGLY(K),K=1,KD)
!     PRINT 704,(SIGLV(K),K=1,KP)
  703 FORMAT(1H ,'PTOP =',F6.2)
  704 FORMAT(1H ,7F10.6)
      DO 913 K=1,KP
       SGLVNU(K) = SIGLV(K)
       IF (K.LE.KD) SIGLNU(K) = SIGLY(K)
  913 CONTINUE
      DO 77 K=1,KD
         Q(K) = SIGLNU(KD+1-K)
   77 CONTINUE
      PSS=    1013250.
      QMH(1)=0.
      QMH(KP)=1.
      DO 1 K=2,KD
      QMH(K)=0.5*(Q(K-1)+Q(K))
1     CONTINUE
      PD(1)=0.
      PD(KP2)=PSS
      DO 2 K=2,KP
      PD(K)=Q(K-1)*PSS
2     CONTINUE
!       call int_get_fresh_handle(retval)
!       close(retval)
!       write(0,*)' before open in CO2O3'
!       open(unit=retval,file=prsmid,form='UNFORMATTED',iostat=ier)
!       write(0,*)' after open1'
!       do k=1,62
!         write(retval)pd(k)
!       enddo
!       close(retval)
      PLM(1)=0.
      DO 3 K=1,KM
      PLM(K+1)=0.5*(PD(K+1)+PD(K+2))
3     CONTINUE
      PLM(KP)=PSS
      DO 4 K=1,KD
      GTEMP(K)=PD(K+1)**0.2*(1.+PD(K+1)/30000.)**0.8/1013250.
4     CONTINUE
      GTEMP(KP)=0.
!+++  WRITE (6,100) (GTEMP(K),K=1,KD)
!+++  WRITE (6,100) (PD(K),K=1,KP2)
!+++  WRITE (6,100) (PLM(K),K=1,KP)
!***TAPES 41,42 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM (PS=1013MB)
!  THE FOLLOWING PUTS P-DATA INTO MB
      DO 11 I=1,KP
      PD(I)=PD(I)*1.0E-3
      PLM(I)=PLM(I)*1.0E-3
11    CONTINUE
      PD(KP2)=PD(KP2)*1.0E-3
!CC         WRITE (41,101) (PD(K),K=1,KP2)
!CC         WRITE (41,101) (PLM(K),K=1,KP)
!CC         WRITE (42,101) (PLM(K),K=1,KP)
      DO 300 K=1,KP2
       T41(K,1) = PD(K)
  300 CONTINUE
      DO 301 K=1,KP
       T41(K,2) = PLM(K)
       T42(K) = PLM(K)
  301 CONTINUE
!***STORE AS PDT,SO THAT RIGHT PD IS RETURNED TO PTZ
      DO 12 I=1,KP2
      PDT(I)=PD(I)
12    CONTINUE
!***SECOND PASS: PSS=810MB,GTEMP NOT COMPUTED
      PSS=0.8*1013250.
      QMH(1)=0.
      QMH(KP)=1.
      DO 201 K=2,KD
      QMH(K)=0.5*(Q(K-1)+Q(K))
201   CONTINUE
      PD(1)=0.
      PD(KP2)=PSS
      DO 202 K=2,KP
      PD(K)=Q(K-1)*PSS
202   CONTINUE
      PLM(1)=0.
      DO 203 K=1,KM
      PLM(K+1)=0.5*(PD(K+1)+PD(K+2))
203   CONTINUE
      PLM(KP)=PSS
!+++  WRITE (6,100) (PD(K),K=1,KP2)
!+++  WRITE (6,100) (PLM(K),K=1,KP)
!***TAPES 43,44 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM(PS=810 MB)
!  THE FOLLOWING PUTS P-DATA INTO MB
      DO 211 I=1,KP
      PD(I)=PD(I)*1.0E-3
      PLM(I)=PLM(I)*1.0E-3
211   CONTINUE
      PD(KP2)=PD(KP2)*1.0E-3
!CC       WRITE (43,101) (PD(K),K=1,KP2)
!CC       WRITE (43,101) (PLM(K),K=1,KP)
!CC       WRITE (44,101) (PLM(K),K=1,KP)
      DO 302 K=1,KP2
       T43(K,1) = PD(K)
  302 CONTINUE
      DO 303 K=1,KP
       T43(K,2) = PLM(K)
       T44(K) = PLM(K)
  303 CONTINUE
!***RESTORE PD
      DO 212 I=1,KP2
      PD(I)=PDT(I)
212   CONTINUE
100   FORMAT (1X,5E20.13)
101   FORMAT (5E16.9)
      RETURN
      END SUBROUTINE SIGP
!---------------------------------------------------------------------
      SUBROUTINE SINTR2
!....
!     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!     REAL P1,P2,PA,TRNSLO,CORE,TRANSA,PATH,UEXP,SEXP,ETA,SEXPV
      COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N
!     COMMON/PRESS/ PA(109)
!     COMMON/TRAN/ TRANSA(109,109)
!     COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP
      DO 70 L=1,109
      IP1=L
      IF (P2-PA(L)) 65,65,70
   70 CONTINUE
   65 I=IP1-1
      IF (IP1.EQ.1) IP1=2
      IF (I.EQ.0) I=1
      DO 80 L=1,109
      JP1=L
      IF (P1-PA(L)) 75,75,80
   80 CONTINUE
   75 J=JP1-1
      IF (JP1.EQ.1) JP1=2
      IF (J.EQ.0) J=1
      JJJ=J
      III=I
      J=JJJ
      JP1=J+1
      I=III
      IP1=I+1
!  DETERMINE ETAP,THE VALUE OF ETA TO USE BY LINEAR INTERPOLATION
!    FOR PETA(=0.5*(P1+P2))
      PETA=P2
      DO 90 L=1,109
      IETAP1=L
      IF (PETA-PA(L)) 85,85,90
90    CONTINUE
85    IETA=IETAP1-1
      IF (IETAP1.EQ.1) IETAP1=2
      IF (IETA.EQ.0) IETA=1
      ETAP=ETA(IETA)+(PETA-PA(IETA))*(ETA(IETAP1)-ETA(IETA))/ &
       (PA(IETAP1)-PA(IETA))
      SEXP=SEXPV(IETA)+(PETA-PA(IETA))*(SEXPV(IETAP1)- &
       SEXPV(IETA))/ (PA(IETAP1)-PA(IETA))
      PIPMPI=PA(IP1)-PA(I)
      UP2P1=(PATH(P2,P1,CORE,ETAP))**UEXP
      IF (I-J) 126,126,127
  126 CONTINUE
      TRIP=(CA(IP1)*DLOG(1.0D0+XA(IP1)*UP2P1))**(SEXP/UEXP)
      TRI=(CA(I)*DLOG(1.0D0+XA(I)*UP2P1))**(SEXP/UEXP)
      TRNSLO=1.0D0-((PA(IP1)-P2)*TRI+(P2-PA(I))*TRIP)/PIPMPI
      GO TO 128
  127 TIJ=TRANSA(I,J)
      TIPJ=TRANSA(I+1,J)
      TIJP=TRANSA(I,J+1)
      TIPJP=TRANSA(I+1,J+1)
      UIJ=(PATH(PA(I),PA(J),CORE,ETAP))**UEXP
      UIPJ=(PATH(PA(I+1),PA(J),CORE,ETAP))**UEXP
      UIJP=(PATH(PA(I),PA(J+1),CORE,ETAP))**UEXP
      UIPJP=(PATH(PA(I+1),PA(J+1),CORE,ETAP))**UEXP
      PRODI=CA(I)*XA(I)
      PRODIP=CA(I+1)*XA(I+1)
      PROD=((PA(I+1)-P2)*PRODI+(P2-PA(I))*PRODIP)/PIPMPI
      XINT=((PA(I+1)-P2)*XA(I)+(P2-PA(I))*XA(I+1))/PIPMPI
      CINT=PROD/XINT
      AIJ=(CINT*DLOG(1.0D0+XINT*UIJ))**(SEXP/UEXP)
      AIJP=(CINT*DLOG(1.0D0+XINT*UIJP))**(SEXP/UEXP)
      AIPJ=(CINT*DLOG(1.0D0+XINT*UIPJ))**(SEXP/UEXP)
      AIPJP=(CINT*DLOG(1.0D0+XINT*UIPJP))**(SEXP/UEXP)
      EIJ=TIJ+AIJ
      EIPJ=TIPJ+AIPJ
      EIJP=TIJP+AIJP
      EIPJP=TIPJP+AIPJP
      DTDJ=(EIJP-EIJ)/(PA(J+1)-PA(J))
      DTDPJ=(EIPJP-EIPJ)/(PA(J+1)-PA(J))
      EPIP1=EIJ+DTDJ*(P1-PA(J))
      EPIPP1=EIPJ+DTDPJ*(P1-PA(J))
      EPP2P1=((PA(I+1)-P2)*EPIP1+(P2-PA(I))*EPIPP1)/PIPMPI
      TRNSLO=EPP2P1-(CINT*DLOG(1.0D0+XINT*UP2P1))**(SEXP/UEXP)
      IF (I.GE.108.OR.J.GE.108) GO TO 350
      IF (I-J-2) 350,350,355
355   CONTINUE
      TIP2J=TRANSA(I+2,J)
      TIP2JP=TRANSA(I+2,J+1)
      TI2J2=TRANSA(I+2,J+2)
      TIJP2=TRANSA(I,J+2)
      TIPJP2=TRANSA(I+1,J+2)
      UIP2J=(PATH(PA(I+2),PA(J),CORE,ETAP))**UEXP
      UIJP2=(PATH(PA(I),PA(J+2),CORE,ETAP))**UEXP
      UIPJP2=(PATH(PA(I+1),PA(J+2),CORE,ETAP))**UEXP
      UI2J2=(PATH(PA(I+2),PA(J+2),CORE,ETAP))**UEXP
      UIP2JP=(PATH(PA(I+2),PA(J+1),CORE,ETAP))**UEXP
      AIJP2=(CINT*DLOG(1.0D0+XINT*UIJP2))**(SEXP/UEXP)
      AIPJP2=(CINT*DLOG(1.0D0+XINT*UIPJP2))**(SEXP/UEXP)
      AIP2J=(CINT*DLOG(1.0D0+XINT*UIP2J))**(SEXP/UEXP)
      AIP2JP=(CINT*DLOG(1.0D0+XINT*UIP2JP))**(SEXP/UEXP)
      AI2J2=(CINT*DLOG(1.0D0+XINT*UI2J2))**(SEXP/UEXP)
      EIP2J=TIP2J+AIP2J
      EIP2JP=TIP2JP+AIP2JP
      EIJP2=TIJP2+AIJP2
      EIPJP2=TIPJP2+AIPJP2
      EI2J2=TI2J2+AI2J2
      CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIJ,EIJP,EIJP2,P1,EI)
      CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIPJ,EIPJP,EIPJP2,P1,EP)
      CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIP2J,EIP2JP,EI2J2,P1,EP2)
      CALL QINTRP(PA(I),PA(I+1),PA(I+2),EI,EP,EP2,P2,EPSIL)
      TRNSLO=EPSIL-(CINT*DLOG(1.0D0+XINT*UP2P1))**(SEXP/UEXP)
  350 CONTINUE
  128 CONTINUE
  205 CONTINUE
      RETURN
      END SUBROUTINE SINTR2
      SUBROUTINE CO2O3(SFULL,SHALF,PPTOP,L,LP1,LP2,RATIN,PSC)
!CCC  PROGRAM CO2O3 = CONSOLIDATION OF A NUMBER OF DAN SCHWARZKOPF,GFDL
!                     CODES TO PRODUCE A FILE OF CO2 HGT DATA
!                     FOR ANY VERTICAL COORDINATE (READ BY SUBROUTINE
!                     CONRAD IN THE GFDL RADIATION CODES)-K.A.C. JUN89.
!NOV89--UPDATED (NOV 89) FOR LATEST GFDL LW RADIATION.....K.A.C.
      LOGICAL                 :: opened
      LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
      CHARACTER*80 errmess
!     integer :: retval,kk,ka,kb
!     character(50) :: co2='co2'
      INTEGER etarad_unit61, etarad_unit62, etarad_unit63,IERROR
      DIMENSION SGTEMP(LP1,2),CO2D1D(L,6),CO2D2D(LP1,LP1,6)
!NOV89
      DIMENSION CO2IQ2(LP1,LP1,6),CO2IQ3(LP1,LP1,6),CO2IQ5(LP1,LP1,6)
!NOV89
      DIMENSION T41(LP2,2),T42(LP1), &
                T43(LP2,2),T44(LP1)
      DIMENSION T20(LP1,LP1,3),T21(LP1,LP1,3)
      DIMENSION T22(LP1,LP1,3),T23(LP1,LP1,3)
      DIMENSION SGLVNU(LP1),SIGLNU(L)
      DIMENSION SFULL(LP1),SHALF(L)
!     DIMENSION STEMP(LP1),GTEMP(LP1)
!     DIMENSION CDTM51(L),CO2M51(L),C2DM51(L)
!     DIMENSION CDTM58(L),CO2M58(L),C2DM58(L)
!     DIMENSION CDT51(LP1,LP1),CO251(LP1,LP1),C2D51(LP1,LP1)
!     DIMENSION CDT58(LP1,LP1),CO258(LP1,LP1),C2D58(LP1,LP1)
!NOV89
!     DIMENSION CDT31(LP1),CO231(LP1),C2D31(LP1)
!     DIMENSION CDT38(LP1),CO238(LP1),C2D38(LP1)
!     DIMENSION CDT71(LP1),CO271(LP1),C2D71(LP1)
!     DIMENSION CDT78(LP1),CO278(LP1),C2D78(LP1)
!     DIMENSION CO211(LP1),CO218(LP1)
!     EQUIVALENCE (CDT31(1),CO2IQ2(1,1,1)),(CO231(1),CO2IQ2(1,1,2))
!     EQUIVALENCE (C2D31(1),CO2IQ2(1,1,3)),(CDT38(1),CO2IQ2(1,1,4))
!     EQUIVALENCE (CO238(1),CO2IQ2(1,1,5)),(C2D38(1),CO2IQ2(1,1,6))
!     EQUIVALENCE (CDT71(1),CO2IQ3(1,1,1)),(CO271(1),CO2IQ3(1,1,2))
!     EQUIVALENCE (C2D71(1),CO2IQ3(1,1,3)),(CDT78(1),CO2IQ3(1,1,4))
!     EQUIVALENCE (CO278(1),CO2IQ3(1,1,5)),(C2D78(1),CO2IQ3(1,1,6))
!     EQUIVALENCE (CO211(1),CO2IQ5(1,1,2)),(CO218(1),CO2IQ5(1,1,5))
!NOV89
!     EQUIVALENCE (STEMP(1),SGTEMP(1,1)),(GTEMP(1),SGTEMP(1,2))
!     EQUIVALENCE (CDTM51(1),CO2D1D(1,1)),(CO2M51(1),CO2D1D(1,2))
!     EQUIVALENCE (C2DM51(1),CO2D1D(1,3)),(CDTM58(1),CO2D1D(1,4))
!     EQUIVALENCE (CO2M58(1),CO2D1D(1,5)),(C2DM58(1),CO2D1D(1,6))
!     EQUIVALENCE (CDT51(1,1),CO2D2D(1,1,1)),(CO251(1,1),CO2D2D(1,1,2))
!     EQUIVALENCE (C2D51(1,1),CO2D2D(1,1,3)),(CDT58(1,1),CO2D2D(1,1,4))
!     EQUIVALENCE (CO258(1,1),CO2D2D(1,1,5)),(C2D58(1,1),CO2D2D(1,1,6))

!
!    Deallocate before reading. This is required for nested domain init.
!
      IF(ALLOCATED (CO251))DEALLOCATE(CO251)
      IF(ALLOCATED (CDT51))DEALLOCATE(CDT51)
      IF(ALLOCATED (C2D51))DEALLOCATE(C2D51)
      IF(ALLOCATED (CO258))DEALLOCATE(CO258)
      IF(ALLOCATED (CDT58))DEALLOCATE(CDT58)
      IF(ALLOCATED (C2D58))DEALLOCATE(C2D58)
      IF(ALLOCATED (STEMP))DEALLOCATE(STEMP)
      IF(ALLOCATED (GTEMP))DEALLOCATE(GTEMP)
      IF(ALLOCATED (CO231))DEALLOCATE(CO231)
      IF(ALLOCATED (CDT31))DEALLOCATE(CDT31)
      IF(ALLOCATED (C2D31))DEALLOCATE(C2D31)
      IF(ALLOCATED (CO238))DEALLOCATE(CO238)
      IF(ALLOCATED (CDT38))DEALLOCATE(CDT38)
      IF(ALLOCATED (C2D38))DEALLOCATE(C2D38)
      IF(ALLOCATED (CO271))DEALLOCATE(CO271)
      IF(ALLOCATED (CDT71))DEALLOCATE(CDT71)
      IF(ALLOCATED (C2D71))DEALLOCATE(C2D71)
      IF(ALLOCATED (CO278))DEALLOCATE(CO278)
      IF(ALLOCATED (CDT78))DEALLOCATE(CDT78)
      IF(ALLOCATED (C2D78))DEALLOCATE(C2D78)
      IF(ALLOCATED (CO2M51))DEALLOCATE(CO2M51)
      IF(ALLOCATED (CDTM51))DEALLOCATE(CDTM51)
      IF(ALLOCATED (C2DM51))DEALLOCATE(C2DM51)
      IF(ALLOCATED (CO2M58))DEALLOCATE(CO2M58)
      IF(ALLOCATED (CDTM58))DEALLOCATE(CDTM58)
      IF(ALLOCATED (C2DM58))DEALLOCATE(C2DM58)
!
      ALLOCATE(CO251(LP1,LP1))
      ALLOCATE(CDT51(LP1,LP1))
      ALLOCATE(C2D51(LP1,LP1))
      ALLOCATE(CO258(LP1,LP1))
      ALLOCATE(CDT58(LP1,LP1))
      ALLOCATE(C2D58(LP1,LP1))
      ALLOCATE(STEMP(LP1))
      ALLOCATE(GTEMP(LP1))
      ALLOCATE(CO231(LP1))
      ALLOCATE(CDT31(LP1))
      ALLOCATE(C2D31(LP1))
      ALLOCATE(CO238(LP1))
      ALLOCATE(CDT38(LP1))
      ALLOCATE(C2D38(LP1))
      ALLOCATE(CO271(LP1))
      ALLOCATE(CDT71(LP1))
      ALLOCATE(C2D71(LP1))
      ALLOCATE(CO278(LP1))
      ALLOCATE(CDT78(LP1))
      ALLOCATE(C2D78(LP1))
      ALLOCATE(CO2M51(L))
      ALLOCATE(CDTM51(L))
      ALLOCATE(C2DM51(L))
      ALLOCATE(CO2M58(L))
      ALLOCATE(CDTM58(L))
      ALLOCATE(C2DM58(L))
      IF ( wrf_dm_on_monitor() ) THEN
        DO i = 61,99
          INQUIRE ( i , OPENED = opened )
          IF ( .NOT. opened ) THEN
            etarad_unit61 = i
            GOTO 2061
          ENDIF
        ENDDO
        etarad_unit61 = -1
 2061   CONTINUE
        DO i = 62,99
          INQUIRE ( i , OPENED = opened )
          IF ( .NOT. opened ) THEN
            etarad_unit62 = i
            GOTO 2062
          ENDIF
        ENDDO
        etarad_unit62 = -1
 2062   CONTINUE
        DO i = 63,99
          INQUIRE ( i , OPENED = opened )
          IF ( .NOT. opened ) THEN
            etarad_unit63 = i
            GOTO 2063
          ENDIF
        ENDDO
        etarad_unit63 = -1
 2063   CONTINUE
      ENDIF
      CALL wrf_dm_bcast_bytes ( etarad_unit61 , IWORDSIZE )
      IF ( etarad_unit61 < 0 ) THEN
        CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' )
      ENDIF
      CALL wrf_dm_bcast_bytes ( etarad_unit62 , IWORDSIZE )
      IF ( etarad_unit62 < 0 ) THEN
        CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' )
      ENDIF
      CALL wrf_dm_bcast_bytes ( etarad_unit63 , IWORDSIZE )
      IF ( etarad_unit63 < 0 ) THEN
        CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' )
      ENDIF
        IF ( wrf_dm_on_monitor() ) THEN
          OPEN(etarad_unit61,FILE='tr49t85',                  &
               FORM='FORMATTED',STATUS='OLD',ERR=9061,IOSTAT=IERROR)
        ENDIF
        IF ( wrf_dm_on_monitor() ) THEN
          OPEN(etarad_unit62,FILE='tr49t67',                  &
               FORM='FORMATTED',STATUS='OLD',ERR=9062,IOSTAT=IERROR)
        ENDIF
        IF ( wrf_dm_on_monitor() ) THEN
          OPEN(etarad_unit63,FILE='tr67t85',                  &
               FORM='FORMATTED',STATUS='OLD',ERR=9063,IOSTAT=IERROR)
        ENDIF

!===>  GET SGTEMP AND OUTPUT WHICH USED TO BE ON UNITS 41,42,43,44....
      LREAD = 0
!     DO KKK=1,L
!JD      READ(23)SIGLNU(KKK)
!      SIGLNU(KKK)=1.-FLOAT(KKK)/LP1
!     END DO
      CALL CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, &
                  SFULL,SHALF,PPTOP,LREAD,L,LP1,LP2)
!--- ETA_PORT: escala de concentracao por pressao (Curtis-Godson): amount x
!    broadening ~ r  <=>  avaliar a tabela 330ppm em pressoes r*p. PSC=ppm/330.
      T41 = T41*PSC
      T42 = T42*PSC
      T43 = T43*PSC
      T44 = T44*PSC
!       call int_get_fresh_handle(retval)
!       close(retval)
!       open(unit=retval,file=co2,form='UNFORMATTED',iostat=ier)
!       do kk=1,2
!         write(retval)(sgtemp(k,kk),k=1,61)
!       enddo
      DO K=1,LP1
        STEMP(K)=SGTEMP(K,1)
        GTEMP(K)=SGTEMP(K,2)
      ENDDO
!===>  INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID..
!         IR=1,IQ=1 IS FOR COMMON /CO2BD3/ IN RADIATION CODE...
!           FOR THE CONSOLIDATED 490-850 CM-1 BAND...
!NOV89
!     ICO2TP=61
      ICO2TP=etarad_unit61
!NOV89
      IR = 1
      RATIO = RATIN
      NMETHD = 2
      CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2)
      IR = 1
      RATIO = RATIN
      NMETHD = 1
      CALL CO2INT(ICO2TP,T41,T42,T20,RATIO,IR,NMETHD,L,LP1,LP2)
      IR = 1
      RATIO = RATIN
      NMETHD = 2
      CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2)
      IR = 1
      RATIO = RATIN
      NMETHD = 1
      CALL CO2INT(ICO2TP,T43,T44,T21,RATIO,IR,NMETHD,L,LP1,LP2)
!===>    FILL UP THE CO2D1D ARRAY
!       THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND
!         THEIR DERIVATIVES FOR TAU(I,I+1),I=1,LEVS,
!         WHERE THE VALUES ARE NOT OBTAINED BY QUADRATURE BUT ARE THE
!         ACTUAL TRANSMISSIVITIES,ETC,BETWEEN A PAIR OF PRESSURES. THESE
!         ARE USED ONLY FOR NEARBY LAYER CALCULATIONS INCLUDING H2O..
!
      IQ = 1
      CALL CO2IN1(T20,T21,CO2D1D,IQ,L,LP1)
!       do kk=1,6
!         write(retval)(co2d1d(k,kk),k=1,60)
!       enddo
      DO K=1,L
        CDTM51(K)=CO2D1D(K,1)
        CO2M51(K)=CO2D1D(K,2)
        C2DM51(K)=CO2D1D(K,3)
        CDTM58(K)=CO2D1D(K,4)
        CO2M58(K)=CO2D1D(K,5)
        C2DM58(K)=CO2D1D(K,6)
      ENDDO
!
!===>    FILL UP THE CO2D2D ARRAY
!    THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND THEIR DERIVATIVES
!        FROM 109-LEVEL LINE-BY-LINE CALCULATIONS MADE USING THE 1982
!        MCCLATCHY TAPE (12511 LINES),CONSOLIDATED,INTERPOLATED
!        TO THE MRF VERTICAL COORDINATE,AND RE-CONSOLIDATED TO A
!        200 CM-1 BANDWIDTH. THE INTERPOLATION METHOD IS DESCRIBED IN
!        SCHWARZKOPF AND FELS (J.G.R.,1985).
!
      CALL CO2INS(T22,T23,CO2D2D,IQ,L,LP1,1)
!       do kk=1,6
!         write(retval)((co2d2d(ka,kb,kk),ka=1,61),kb=1,61)
!       enddo
      DO K1=1,LP1
      DO K2=1,LP1
        CDT51(K1,K2)=CO2D2D(K1,K2,1)
        CO251(K1,K2)=CO2D2D(K1,K2,2)
        C2D51(K1,K2)=CO2D2D(K1,K2,3)
        CDT58(K1,K2)=CO2D2D(K1,K2,4)
        CO258(K1,K2)=CO2D2D(K1,K2,5)
        C2D58(K1,K2)=CO2D2D(K1,K2,6)
      ENDDO
      ENDDO
!
!NOV89
!===>  INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID..
!         IR=2,IQ=2 IS FOR COMMON /CO2BD2/ IN RADIATION CODE...
!           FOR THE CONSOLIDATED 490-670 CM-1 BAND...
!     ICO2TP=62
      ICO2TP=etarad_unit62
      IR = 2
      RATIO = RATIN
      NMETHD = 2
      CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2)
      CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2)
      IQ = 2
      CALL CO2INS(T22,T23,CO2IQ2,IQ,L,LP1,2)
!       do kk=1,6
!         write(retval)(co2iq2(k,1,kk),k=1,61)
!       enddo
      DO K=1,LP1
        CDT31(K)=CO2IQ2(K,1,1)
        CO231(K)=CO2IQ2(K,1,2)
        C2D31(K)=CO2IQ2(K,1,3)
        CDT38(K)=CO2IQ2(K,1,4)
        CO238(K)=CO2IQ2(K,1,5)
        C2D38(K)=CO2IQ2(K,1,6)
      ENDDO
!===>  INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID..
!         IR=3,IQ=3 IS FOR COMMON /CO2BD4/ IN RADIATION CODE...
!           FOR THE CONSOLIDATED 670-850 CM-1 BAND...
!     ICO2TP=63
      ICO2TP=etarad_unit63
      IR = 3
      RATIO = RATIN
      NMETHD = 2
      CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2)
      CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2)
      IQ = 3
      CALL CO2INS(T22,T23,CO2IQ3,IQ,L,LP1,3)
!       do kk=1,6
!         write(retval)(co2iq3(k,1,kk),k=1,61)
!       enddo
!       close(retval)
      DO K=1,LP1
        CDT71(K)=CO2IQ3(K,1,1)
        CO271(K)=CO2IQ3(K,1,2)
        C2D71(K)=CO2IQ3(K,1,3)
        CDT78(K)=CO2IQ3(K,1,4)
        CO278(K)=CO2IQ3(K,1,5)
        C2D78(K)=CO2IQ3(K,1,6)
      ENDDO
!---      FOLLOWING CODE NOT WORKING AND NOT NEEDED YET
!===>  INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID..
!         IR=4,IQ=5 IS FOR COMMON /CO2BD5/ IN RADIATION CODE...
!           FOR THE 4.3 MICRON BAND...
! NOT USED YET      ICO2TP=65
! NOT USED YET      IR = 4
! NOT USED YET      RATIO = 1.0
! DAN SCHWARZ --- USE 300PPMV  RATIO = 0.9091   (NOT TESTED YET).....
! NOT USED YET      NMETHD = 2
! NOT USED YET      CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD)
! NOT USED YET      CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD)
! NOT USED YET      IQ = 5
! NOT USED YET      CALL CO2INS(T22,T23,CO2IQ5,IQ)
!NOV89
!...     WRITE DATA TO DISK..
!            ...SINCE THESE CODES ARE COMPILED WITH AUTODBL,THE CO2 DATA
!               IS CONVERTED TO SINGLE PRECISION IN A LATER JOB STEP..

! NOT USED YET      WRITE(66) CO211
! NOT USED YET      WRITE(66) CO218
!NOV89
         IF ( wrf_dm_on_monitor() ) THEN
           CLOSE (etarad_unit61)
           CLOSE (etarad_unit62)
           CLOSE (etarad_unit63)
         ENDIF

      RETURN
9061 CONTINUE
     WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr49t85 on unit ',etarad_unit61
     write(0,*)' IERROR=',IERROR
     CALL wrf_error_fatal(errmess)
9062 CONTINUE
     WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr49t67 on unit ',etarad_unit62
     write(0,*)' IERROR=',IERROR
     CALL wrf_error_fatal(errmess)
9063 CONTINUE
     WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr67t85 on unit ',etarad_unit63
     write(0,*)' IERROR=',IERROR
     CALL wrf_error_fatal(errmess)
      END SUBROUTINE CO2O3
