C=======================================================================
C=======================================================================
C
C        GODDARD SPACE FLIGHT CENTER
C
C        AIRS
C
C        UMBC Rapid Transmittance Algorithm (RTA) Wrapper 
C
!F77====================================================================
C
!ROUTINE NAME: INIT_RTA
C
!ABSTRACT:
C
C    Calculate the fast transmittance code predictors for a given profile.
C    The predictors depend on gas amount, temperature, and secant angle.
C    Temperature ratios and moisture moisture amounts are converted to
C    offsets relative to a reference profile are calculated.
C
!CALL PROTOCOL:
C
C    SUBROUTINE init_rta( mode, iphys, NLAY, Psurf, Tsurf,
C                          Temp, wcd, ocd, wlcd, Pobs, Tlayer) 
C
C
!INPUT PARAMETERS:
C
C   Name        Type     DESCRIPTION
C   ----        ----     -----------
C   mode      integer    bit-wise parameter to direct computations
C   iphys     integer    backward compatibility paramerer (GSFC version only)
C   NLAY      integer    # of layers in this profile.
C   Psurf     real       surface pressure (milli-bar)
C   Tsurf     real       surface temperature (K)
C   Temp(L)   real       temperature at level == Pobs(L)
C   wcd(L)    real       layer column density of water
C                        between Pobs(L) and Pobs(L-1)
C                        units are molecules/cm^2
C                        top level Pobs(0) assumed equal to 0.005
C   ocd(L)    real       layer column density of ozone
C   wlcd(L)   real       layer column density of liquid water (NOT used for IR)
C   Pobs(L)   real       pressure at levels, milli-bar
C
!OUTPUT PARAMETERS:
C
C   Name        Type     DESCRIPTION
C   ----        ----     -----------
C   Tlayer      real     layer temperature (to be passed to AIRSRAD)
C
C   all other output is sent to common block RTA_PROF.COM
C
!INPUT/OUTPUT PARAMETERS:  None
C
!RETURN VALUES: none
C
!PARENT(S):
C
C    any routine computing transmittances
C
!ROUTINES CALLED:
C 
C    softexit('INIT_RTA',) if problems encountered
C
!FILES ACCESSED:  None
C
!COMMON BLOCKS:
C
C    RTA_STAT.COM
C    RTA_PROF.COM
C
!DESCRIPTION:
C
C
C
!ALGORITHM REFERENCES: 
C 
C    1) AIRS ATBD
C    2) UMBC routines calpar.f and sunpar.f  (5/98 version)
C    3) Scott Hannon, L. Larrabee Strow, W. Wallace McMillan
C       "Atmospheric Infrared Fast Transmittance Models: A Comparison
C        of Two Approaches." 1996  S.P.I.E. vol 2830, p.94 
C
!KNOWN BUGS AND LIMITATIONS:  None
C
!ROUTINE HISTORY:
C   Date      Programmer      Description
C   Aug. 98   Scott Hannon    Original Package Delivery as CALPAR, SUNPAR
C   Sep. 98   Chris Barnet    merged code to form INIT_RTA()
C   Nov. 98   C. Barnet       added Tlayer to call list
C   12/98     C.Barnet        make solzang test .ge. solzang_limit (was 90.0)
C
!END====================================================================
C
C
C=======================================================================
c
c  Documentation to summarize models supported by this library.
c
c  Instruments supported at this time are:
c    AIRS PDR RTA
c    CrIS 3/98 L=0.8 version
c    CrIS 10/98 L=0.8, 1.0, 1.2 versions
c    IASI 8/98 version
c    NAST 9/98 version
c
c
c  -------------------------------------------------------
c                      AIRS 3/98 RTA
c  -------------------------------------------------------
c
c             AIRS  C  F  W  O CO CH4 TOT N*TOT
c   1 FWO     1561  4  8 11  5  0   0  28 43708
c   2 FOW      180  4  8 11 10  0   0  33  5940
c   3 FMW      122  4  8 11  0  0   9  32  3904
c   4 FCOW     108  4 11 13  3 11   0  42  4536
c   5 FWO-bfsw 131  4 11  3  1  0   0  19  2489 
c   6 FWO-mfmw 144  4  8  7  1  0   0  20  2880
c   7 FWO-mfbw 127  4  8 13  1  0   0  26  3302
c             ----                        -----
c             2373                        66759
c
c    bfsm ==  big   fixed,  small water
c    mfmm == medium fixed, medium water
c    mfbm == medium fixed,   big  water
c
c   66759 * 100 levels, real*4 -->  26,703,600 bytes
c
cc OPTRAN coefficients
c  419 channels, 300 layers, 9 COEFS, 4 PROFILE AVERAGES?
c  INDH2O  2373
c  WAZOP    419
c  WAVGOP  4*300            ==> 4,541,168b
c  COFH2O  9*300*419
c
c CO2 coefficients
c  INDCO2 = 2373            ==>     9,492b
c  COFCO2 =  954*4*100      ==> 1,526,400b
c
c DOWNWELLING: COEFF, LABOVE
c  5*2373 = 11865 + 2373 = 14238 ==> 56,952b
c
c
c ----------------------------------------------------------------------
c                           CrIS 1/98 RTA
c ----------------------------------------------------------------------
c
c  - L=0.79994, 761 chl's from  630.05 (#1008) to 1105.08 (#1778)
c  - L=0.39997  481 chl's from 1205.09 ( #964) to 1805.13 (#1444)
c  - L=0.20008  161 chl's from 2104.15 ( #842) to 2503.99 (#1002)
c
c    CrIS RTA
c
c          CrIS  C  F  W  O CO CH4 TOT N*TOT
c   8 FWO  1105  4  8 11  5  0   0  28 19890 (uses #1's coef memory)
c   9 FOW   137  4  8 11 10  0   0  33  4521 (uses #2's coef memory)
c  10        44  4  8 11 11  0   0  34  1496 (uses #4's coef memory)
c  11 FWO-b 117  4 11  3  1  0   0  19  2223 (uses #5's coef memory)
c          ----                        -----
c          1403                        28130
c
c 100 levels, real*4 --> 11,252,000 bytes
c
c  
c  -------------------------------------------------------
c                   AIRS 12/95 RTA
c  -------------------------------------------------------
c
c   12/95  RTA
c          AIRS  C  F  W  O CO CH4 TOT N*TOT
c   0      2380  4  8 13  9  0   0  34 80920
c
c 2380*34 = 80920 * 100 levels, real*4 --> 32,368,000
c  108*42    4536 * 100 levels, real*4 -->  1,814,400
c                                          ----------
c  (3/98 is 26,703,600 bytes)              34,182,400 bytes
c
c 2380*34 = 80920 * 100 levels, real*4 --> 32,368,000
c  150*42 =  6300 * 100 levels, real*4 -->  2,520,000
c                                          ----------
c                                          34,888,000 bytes
c 
c
c 2373*42 = 99666 * 100 levels, real*4 --> 39,866,400 bytes
c
c ----------------------------------------------------------------------
c
c  ------------------------------------------------------------------------
c  ------------------------------------------------------------------------
c                  10/98 CrIS RTA at L=0.8, 1.0, 1.2
c                              IASI RTA 
c                              NAST RTA 
c  ------------------------------------------------------------------------
c
c  NAST
c  ----
c
c   from NAST_FTC_Aug98/Data/Lists/list_band?.txt.  The original
c   UMBC delivery (8/98) had 4149 channels in band.2, of which, 853
c   between 1974.39 and 2180.03 were deleted because they were duplicates.
c   All the deleted channels were from the original set #3 (now calt14)
c   which had 6172 channels and now has 5319 channels.  The original
c   set had 11305.  The program cnv_iasi.F was used to delete the channels.
c
c   band #1 has 3526 channels    605.095 to 1454.879  0.241074  2.07406
c   band #2 has 3296 channels   1180.055 to 1974.393  0.241074  2.07406
c   band #3 has 3630 channels   1955.107 to 2829.963  0.241074  2.07406
c              -----
c              10452
c
c           65 layers, real*4 --> 86,204,300 bytes
c
c     NAST    100
c              35  p(35) = 51.528
c       1 <->  36  p(36) = 56.126
c      62 <->  97 
c      65 <-> 100
c
c
c  IASI
c  ----
c    605.095-2829.96   9230 channels  L = 2.0   Gaussian Apodized
c
c    f = 605.095 + dv*(n-1),   dv=0.241074,  L = 2.07406
c
c  CrIS
c  ----
c    
c
c               band.1 band.2 band.3  total
c    CrIS L=0.8    939    483    294   1716
c    CrIS L=1.0   1171    601    367   2139
c    CrIS L=1.2   1408    723    441   2572
c
c  - L=0.79994  761 chl's from  630.05 (#1008) to 1105.08 (#1778)
c  - L=0.39997  481 chl's from 1205.09 ( #964) to 1805.13 (#1444)
c  - L=0.20008  161 chl's from 2104.15 ( #842) to 2503.99 (#1002)
c
c
c  CrIS L=0.8
c - L=0.80210  939 chls from   605.29(# 971) to   1190.00(#1909) dv=0.6233637
c - L=0.40100  483 chls from  1199.50(# 962) to   1800.50(#1444) dv=1.2468829
c - L=0.20050  294 chls from  2099.75(# 842) to   2830.42(#1135) dv=2.4937649
c
c  CrIS L=1.0
c - L=1.00060 1171 chls from   605.14(#1211) to   1189.79(#2381) 0.4997002
c - L=0.50030  601 chls from  1200.28(#1201) to   1799.92(#1801) 0.9994004
c - L=0.25010  367 chls from  2099.16(#1050) to   2830.87(#1416) 1.9992002
c
c  CrIS L=1.2
c - L=1.20310 1408 chls from   605.10(#1456) to   1189.84(#2863) 0.4155930
c - L=0.60160  723 chls from  1200.13(#1444) to   1800.20(#2166) 0.8311171
c - L=0.30080  441 chls from  2099.40(#1263) to   2830.78(#1703) 1.6622342
c
c
c
c        new   orig                       # of channels
c      calt#  calt#  C  F  W  O CO tot  L=0.8  L=1.0  L=1.2   IASI   NAST
c         12      1  4 11  7  5  -  27    729    940   1166   4371   3730
c         13      2  4  8 11 10  -  33    236    301    358    737    621
c         14      3  4  8 13  3  -  28    702    837    974   3405   5319
c         15      4  4 11 13  9 11  48     49     61     74    717    782
c                                       -----  -----  -----  -----  -----
c                                        1716   2139   2572   9230  10452
c
c  ------------------------------------------------------------------------
c  ------------------------------------------------------------------------
c  ------------------------------------------------------------------------
c
c                         table to go from FPREDX(i,L) --> FPRED(j,L)
c      
c                               values of i are given in table
c
c                       ........... AIRS .......... | ..... CrIS .... |  .... NAST ....
c        j   predictor  X=1   2   3   4   5   6   7 |   8   9  10  11 |  12  13  14  15
c      ----- --------- ---- --- --- --- --- --- --- | --- --- --- --- | --- --- --- ---
c        1   a            1   1   1   1   1   1   1 |   1   1   1   1 |   1   1   1   1
c        2   a^2          2   2   2   2   2   2   2 |   4   4   4   4 |   4   4   4   4
c        3   a*Tr         3   3   3   3   3   3   3 |   3   3   3   3 |   3   3   3   3
c        4   a*Tr^2       4   4   4   4   4   4   4 |   5   5   5   5 |   5   5   5   5
c        5   Tr           5   5   5   5   5   5   5 |   2   2   2   2 |   2   2   2   2
c        6   Tr^2         6   6   6   6   6   6   6 |   6   6   6   6 |   6   6   6   6
c        7   a*Trz        7   7   7   7   7   7   7 |   7   7   7   7 |   7   7   7   7
c        8   a*Trz/Tr     8   8   8   .   8   .   . |   8   8   .   . |   8   8   .   8
c        9   a^2*Tr       .   .   .   9   9   .   . |   .   .   .   9 |   9   .   .   9
c       10   a^3          .   .   .  10   .   .   . |   .   .   .   . |   .   .   .   .
c       11   sqrt(a)      .   .   .  11  10   8   8 |   .   .   8  10 |  10   .   8  10
c       12   a^2*Trz      .   .   .   8   .   .   . |   .   .   .   . |   .   .   .   .
c       13   Trz          .   .   .   .  11   .   . |   .   .   .  11 |  11   .   .  11
c       14   sqrt(Tr)     .   .   .   .   .   .   . |   .   .   .   8 |   .   .   .   .
c       15                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       16                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .

c      ======================================================================
c      
c                         table to go from CPRED4(i,L) --> CPRED(j,L)
c      
c                               values of i are given in table
c      
c                       ........... AIRS .......... | ..... CrIS .... |  .... NAST ....
c        j   predictor  X=1   2   3   4   5   6   7 |   8   9  10  11 |  12  13  14  15
c      ----- --------- ---- --- --- --- --- --- --- | --- --- --- --- | --- --- --- ---
c        1   C*a          .   .   .   1   .   .   . |   .   .   .   . |   .   .   .   1
c        2   sqrt(C*a)    .   .   .   2   .   .   . |   .   .   .   . |   .   .   .   2
c        3   C*a*dT       .   .   .   3   .   .   . |   .   .   .   . |   .   .   .   3
c        4   (C*a)^2      .   .   .   4   .   .   . |   .   .   .   . |   .   .   .   4
c        5   C*a*C/Cz     .   .   .   5   .   .   . |   .   .   .   . |   .   .   .   5
c        6   sqrt(C*a)*dT .   .   .   6   .   .   . |   .   .   .   . |   .   .   .   6
c        7   root^4(C*a)  .   .   .   7   .   .   . |   .   .   .   . |   .   .   .   7
c        8 sqrt(C*a)*C/Cz .   .   .   8   .   .   . |   .   .   .   . |   .   .   .   8
c        9   C            .   .   .   9   .   .   . |   .   .   .   . |   .   .   .   9
c       10   a^2*C        .   .   .  10   .   .   . |   .   .   .   . |   .   .   .  10
c       11   a*sqrt(C*a)  .   .   .  11   .   .   . |   .   .   .   . |   .   .   .  11
c       12   Cz           .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       13                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       14                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       15                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       16                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c      
c      ======================================================================
c      
c                         table to go from MPRED3(i,L) --> MPRED(j,L)
c      
c                               values of i are given in table
c      
c                       ........... AIRS .......... | ..... CrIS .... |  .... NAST ....
c        j   predictor  X=1   2   3   4   5   6   7 |   8   9  10  11 |  12  13  14  15
c      ----- --------- ---- --- --- --- --- --- --- | --- --- --- --- | --- --- --- ---
c        1   M*a          .   .   1   .   .   .   . |   .   .   .   . |   .   .   .   .
c        2   sqrt(M*a)    .   .   2   .   .   .   . |   .   .   .   . |   .   .   .   .
c        3   M*a*dT       .   .   3   .   .   .   . |   .   .   .   . |   .   .   .   .
c        4   (M*a)^2      .   .   4   .   .   .   . |   .   .   .   . |   .   .   .   .
c        5   M*a^2        .   .   5   .   .   .   . |   .   .   .   . |   .   .   .   .
c        6   Mz*a         .   .   6   .   .   .   . |   .   .   .   . |   .   .   .   .
c        7   M*dT         .   .   7   .   .   .   . |   .   .   .   . |   .   .   .   .
c        8   TMz*a        .   .   8   .   .   .   . |   .   .   .   . |   .   .   .   .
c        9   sqrt(Mz*a)   .   .   9   .   .   .   . |   .   .   .   . |   .   .   .   .
c       10   M            .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       11   Mz           .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       12                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       13                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       14                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       15                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       16                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c      
c      
c      ======================================================================
c      
c      
c      
c                         table to go from OPRED3(i,L) --> OPRED(j,L)
c      
c                               values of i are given in table
c      
c                       ........... AIRS .......... | ..... CrIS .... |  .... NAST ....
c        j   predictor  X=1   2   3   4   5   6   7 |   8   9  10  11 |  12  13  14  15
c      ----- --------- ---- --- --- --- --- --- --- | --- --- --- --- | --- --- --- ---
c        1   O*a          1   1   .   1   1   1   1 |   1   1   1   1 |   1   1   1   1
c        2   sqrt(O*a)    2   2   .   2   .   .   . |   2   2   2   . |   2   2   2   2
c        3   O*a*dT       3   3   .   3   .   .   . |   3   3   3   . |   3   3   3   3
c        4   (O*a)^2      4   4   .   .   .   .   . |   4   4   4   . |   4   4   .   4
c        5   sqrt(O*a)*dT 5   5   .   .   .   .   . |   5   5   5   . |   .   5   .   5
c        6   O*a*O/Ox     .   6   .   .   .   .   . |   .   6   6   . |   .   6   .   6
c        7 sqrt(O*a)*O/Ox .   7   .   .   .   .   . |   .   7   7   . |   .   7   .   .
c        8   O*a*Oz/Ox    .   8   .   .   .   .   . |   .   8   8   . |   .   .   .   .
c        9 O*a*sqrt(Ox*a) .   9   .   .   .   .   . |   .   9   9   . |   .   9   .   .
c       10   O*a*TOz*a    .  10   .   .   .   .   . |   .  10  10   . |   .  10   .   .
c       11  (O*a)^3/2     .   .   .   .   .   .   . |   .   .  11   . |   .   .   .   9
c       12  a^2*O         .   .   .   .   .   .   . |   .   .   .   . |   5   .   .   .
c       13  (O*a)^2*O/Ox  .   .   .   .   .   .   . |   .   .   .   . |   .   8   .   .
c       14  (O*a)^3       .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   7
c       15  O             .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   8
c       16                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c
c      ======================================================================
c      
c      
c                         table to go from WPRED3(i,L) --> WPRED(j,L)
c      
c                               values of i are given in table
c      
c                       ........... AIRS .......... | ..... CrIS .... |  .... NAST ....
c        j   predictor  X=1   2   3   4   5   6   7 |   8   9  10  11 |  12  13  14  15
c      ----- --------- ---- --- --- --- --- --- --- | --- --- --- --- | --- --- --- ---
c        1   W*a          1   1   1   1   1   1   1 |   1   1   1   1 |   1   1   1   1
c        2   sqrt(W*a)    2   2   2   3   .   .   . |   2   .   .   . |   .   .   2   3
c        3   W*a*W/Wz     3   8   3   8   .   .   8 |   3   8   8   . |   .   8   8   8
c        4   W*a*dT       4   3   4   4   3   3   3 |   4   3   3   3 |   3   3   3   4
c        5   (W*a)^2      5   5   5   5   .   4   4 |   5   5   5   . |   4   5   4   5
c        6  sqrt(W*a)*dT  6   7   6   6   .   .   . |   6   .   .   . |   .   .   5   6
c        7  root^4(W*a)   7   6   7   7   .   .   . |   7   .   .   . |   .   .  10   7
c        8 sqrt(W*a)*W/Wz 8  11  10  12   .   .   . |   8   .   .   . |   .   .   9  12
c        9   (W*a)^3      9   9   8  10   .   6   6 |   9   9   9   . |   .   9   6  10
c       10   W           10   .   9   2   .   .   . |  10  10  10   . |   .   .   .   2
c       11  W*a*dT*|dT|  11   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       12  W*a*Ox*a      .   4   .   .   .   .   . |   .   4   4   . |   .   4   .   .
c       13  W*a*(Ox*a)^2  .  10   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       14  W*a*Mz*a      .   .  11   .   .   .   . |   .   .   .   . |   .   .   .   .
c       15  W*a^2         .   .   .   9   .   7   7 |   .   .   .   . |   7  10   7   9
c       16  W*a*Cz*a      .   .   .  11   .   .   . |   .   .   .   . |   .   .   .  11
c       17  W*a^2*dT      .   .   .  13   .   .   . |   .   .   .   . |   .   .   .  13
c       18 (W*a)^3/2      .   .   .   .   2   2   2 |  11   2   2   2 |   2   2   .   .
c       19 (W*a)^3/2*dT   .   .   .   .   .   5   5 |   .   7   7   . |   5   7   .   .
c       20 (W*a)^3/2*W/Wz .   .   .   .   .   .   9 |   .  11  11   . |   .  11   .   .
c       21  (W*a)^5/4     .   .   .   .   .   .  10 |   .   6   6   . |   6   6   .   .
c       22  (W*a)^2*W/Wz  .   .   .   .   .   .  11 |   .   .   .   . |   .   .  11   .
c       23  W^2*a         .   .   .   .   .   .  12 |   .   .   .   . |   .   .  13   .
c       24 (W*a)^7/4      .   .   .   .   .   .  13 |   .   .   .   . |   .   .   .   .
c       25 (W*a)^3/4      .   .   .   .   .   .   . |   .   .   .   . |   .   .  12   .
c       26                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       27                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       28                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       29                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c       30                .   .   .   .   .   .   . |   .   .   .   . |   .   .   .   .
c      
c      ======================================================================
c      
C    -----
C    where:
C    "a" is the secant of the viewing angle SECANG
C    "Tr" is the temperature ratio PTEMP/RTEMP
C    "Trz" is the pressure weighted temperature ratio above, i.e.
C      the sum i=2 to i=L of { P(i) * ( P(i) -  P(i-1) )* Tr(i-1) }
C      where "P" is the pressure PRES and "L" is the layer number, and
C      Trz(L=1)=0
C    "W" is the water amount ratio PWAMNT/RWAMNT
C    "dT" is the temperature offset PTEMP-RTEMP
C    "Wz" is the pressure weighted water amount above ratio, the
C      sum i=1 to i=L of { P(i) * ( (P(i)-P(i-1) ) * PWAMNT(i) },
C      divided by the same sum except using RWAMNT instead of PWAMNT.
C      For these sums, term P(0) is defined as P(0)=2*P(1) - P(2).
C    "O" is the ozone amount ratio POAMNT/ROAMNT
C    "Oz" is the pressure weighted ozone amount above ratio, the
C      sum i=1 to i=L of { P(i) * ( (P(i)-P(i-1) ) * POAMNT(i) },
C      divided by the same sum except using ROAMNT instead of POAMNT.
C      For these sums, term P(0) is defined as P(0)=2*P(1) - P(2)
C    "Ox" is the unweighted ozone amount above ratio, the
C      sum i=1 to i=L of { POAMNT(i) },
C      divided by the same sum except using ROAMNT instead of POAMNT.
C      For these sums, term P(0) is defined as P(0)=2*P(1) - P(2).
C    "TOz" is the pressure and ozone weighted temperature ratio above,
C      sum i=2 to i=L of { P(i) * ( P(i)-P(i-1) )* dT(i-1) * O(i-1) }
C      and TOz(L=1)=0
C    "C" is the carbon monoxide amount ratio POAMNT/ROAMNT
C    "Cz" is the pressure weighted CO amount above ratio, the
C      sum i=1 to i=L of { P(i) * ( (P(i)-P(i-1) ) * PCAMNT(i) },
C      divided by the same sum except using RCAMNT instead of PCAMNT.
C      For these sums, term P(0) is defined as P(0)=2*P(1) - P(2).
C    "M" is the methane amount ratio PMAMNT/RMAMNT
C    "Mz" is the pressure weighted methane amount above ratio, the
C      sum i=1 to i=L of { P(i) * ( (P(i)-P(i-1) ) * PMAMNT(i) },
C      divided by the same sum except using RMAMNT instead of PMAMNT.
C      For these sums, term P(0) is defined as P(0)=2*P(1) - P(2).
C    "TMz" is the pressure and methane weighted temperature ratio above,
C      sum i=2 to i=L of { P(i) * ( P(i)-P(i-1) )* Tr(i-1) * M(i-1) }
C      and TMz(L=1)=0
C
C    ===================================================================

       subroutine init_rta( mode, iphys, NLAY, Psurf, Tsurf,
     +   Temp, wcd, ocd, wlcd, Pobs, Tlayer ) 
       IMPLICIT NONE


#include "rta_stat.com"
#include "rta_prof.com"

      integer*4   solzang_limit
      parameter   (solzang_limit = 85.0)

C    input variables:
C    ----------------
       integer*4  mode     ! which predictors to calculate
       integer*4  iphys    ! revision level
       integer*4  NLAY     ! # of levels
       real*4     Psurf    ! surface pressure
       real*4     Tsurf    ! surface temperature  
       real*4     Temp(*)  ! temperature profile  > TEMP
       real*4     wcd(*)   ! water profile        > WAMNT
       real*4     ocd(*)   ! ozone profile        > OAMNT
       real*4     wlcd(*)  ! liquid water profile > NOT USED 
       real*4     Pobs(*)  ! pressure for Temp, wcd, ocd, wlcd

C    output variables:
C    -----------------
      real*4      Tlayer(*) ! layer temperature (to be passed to AIRSRAD)
C
C    Local Variables
C    ---------------
       INTEGER*4  I
       INTEGER*4  L
       real*4     a, a2   ! local secant angle, SECANG(L) or SECSUN(L)
       REAL*4    PDP
       REAL*4  PNORM
       REAL*4     DT
       REAL*4     TR
       REAL*4     TZ, TRZ
       REAL*4    A_F
       REAL*4    A_W
       REAL*4  WZREF
       REAL*4     WZ
       REAL*4   AZ_W
       REAL*4    A_O
       REAL*4  XZREF
       REAL*4     XZ
       REAL*4   XZ_O
       REAL*4  OZREF
       REAL*4     OZ
       REAL*4   AZ_O
       REAL*4    TOZ
       REAL*4  TAZ_O
       REAL*4    A_C
       REAL*4     CZ
       REAL*4  CZREF
       REAL*4   AZ_C
       REAL*4    A_M
       REAL*4  MZREF
       REAL*4     MZ
       REAL*4   AZ_M
       REAL*4    TMZ
       REAL*4  TAZ_M
       REAL*4 TJUNKS
       REAL*4 WJUNKA, WJUNKR, WJUNKS, WJUNKZ, WJUNK4
       REAL*4 OJUNKA, OJUNKR,         OJUNKZ,        OJUNKX
       REAL*4 CJUNKA, CJUNKR, CJUNKS, CJUNKZ
       REAL*4 MJUNKA, MJUNKR,         MJUNKZ
C-- OPTRAN:
       INTEGER*4 LL  
       INTEGER*4 LOP
       INTEGER*4 LOPL
       INTEGER*4 LOPU
       INTEGER*4 LU
       REAL*4    WAZ(MAXLAY)
       REAL*4    WAZSUM
       REAL*4    WPZSUM
       REAL*4    WTZSUM
       REAL*4    PZ(MAXLAY)
       REAL*4    DA
       REAL*4    TZopt(MAXLAY)
       REAL*4    POP
       REAL*4    TOP
       REAL*4    PZOP
       REAL*4    TZOP
       REAL*4    ANGOP
       LOGICAL*4   LAST

       real*4   kmole
       parameter (kmole = 1.0/6.0221367E+26)

      logical*4  exactP
      real*4     x
      real*4     Tavg(MAXLAY)

c     variables for the fixed gas multiplied (FIXMUL)
c     -----------------------------------------------
C     FIXMUL: the not-quite-fixed "fixed" amount multiplier:
C     FIXMUL = PFAMNT/RFAMNT  for UMBC profiles where FIXED GAS amount is know
C     This value should be close to (within a few percent of) unity.
C     For the input to init_rta() the water partial pressure is computed
c     as follows:
c     wcd(L) = layer column density in molecules/cm^2 = Na*1000*kilo-moles/cm^2
c     ALT(L) = altitude used in secant computation (kilo-meters)
c     N = wcd(L)/(1.0E5*ALT(L)) == molecules/cm^3
c
c     Na = 6.022045E+23 molecules/mole
c     Pwater = (N/Na)*R*T/Pstd = partial pressure in atm
c
c     N0 = Loschmidt's number = 2.6867E+19 molecules/cm^3
c     N0 = Pstd*Na/R/Tstd,  Pstd=101325E6 dynes/cm^2  Tstd=273.15
c     R/(Na*Pstd) = 1/(N0*Tstd)
c     Pwater = N*T/(N0*Tstd),  in atm

      real*4 dz          ! thickness of layer in cm's
      real*4 Pwater      ! partial pressure of water (atm)
      real*4 Pmult       ! correction factor (UMBC)
      parameter (pmult=0.622)
      real*4 fx(MAXLAY)  ! correction factor (UMBC)
      data fx/
     &1.0031478, 1.0030850, 1.0034791, 1.0032465, 1.0027664, 1.0030644,
     &1.0027682, 1.0024457, 1.0024817, 1.0019743, 1.0021864, 1.0019980,
     &1.0016734, 1.0016632, 1.0015243, 1.0014169, 1.0011749, 1.0012166,
     &1.0007437, 1.0010513, 1.0007979, 1.0008730, 1.0004258, 1.0006266,
     &1.0007262, 1.0003718, 1.0002833, 1.0004674, 1.0004928, 1.0003764,
     &1.0001269, 1.0003734, .99998352, 1.0006066, 1.0003759, 1.0006079,
     &.99989271, .99955584, .99950482, .99980806, 1.0018658, 1.0027118,
     &1.0007736, 1.0016782, 1.0016151, 1.0003560, .99985240, .99975007,
     &.99996031, 1.0001826, 1.0002389, .99999216, 1.0001728, .99990960,
     &1.0003671, .99986994, 1.0000821, 1.0002300, 1.0000069, 1.0001705,
     &1.0000816, .99999931, 1.0002539, 1.0004827, 1.0001965, 1.0003765,
     &1.0000926, 1.0002179, 1.0001384, 1.0001516, 1.0001784, 1.0001839,
     &1.0002541, 1.0003887, 1.0005754, 1.0001825, .99973361, .99969221,
     &1.0005867, 1.0025437, 1.0020681, 1.0005591, 1.0007326, 1.0004324,
     &1.0003146, 1.0006219, 1.0005623, 1.0003755, 1.0008667, 1.0001547,
     &1.0007129, .99958034, .99690980, .99896196, .99937805, .99848837,
     &.99833538, .99744689, .99633090, .99218102/

C     Assign the bottom layer and average temperature
C     -----------------------------------------------

      LBOT = NLAY  ! (1-100)
      BLMULT = (Psurf-Pobs(NLAY-1))/(Pobs(NLAY)-Pobs(NLAY-1))

C     layer tranmittance will be adjusted by BLMULT at end, extrapolated
C     temperature is used for now (..... will be changed before delivery)
      Tavg(LTOP) = Temp(LTOP)
      do L = LTOP+1, NLAY
        Tavg(L) = 0.5*(Temp(L)+Temp(L-1))
      enddo

      call meantemp(NLAY,Pobs,Psurf,Temp, Tlayer)

C     Check stuff if first time to make sure all is OK
C     ------------------------------------------------
      if( firstcall ) then  !firstcall of every profile, set=T by setang
         
        if( iphys.lt.4 ) then
          print 100, iphys
          i=11
          call softexit('INIT_RTA',i)
        endif

        if( NLAY.gt.MAXLAY ) then  ! or ntrlev_ir
          print 400, NLAY, ntrlev_ir, Psurf, Pobs(ntrlev_ir)
          i=8
          call softexit('INIT_RTA',i)
        endif

        exactP = .true.
        do L = LTOP, NLAY
          x = ABS( Pobs(L)-Pstd(L) )
          if( x.gt.2.0E-3 ) then    ! ASCII is f10.3
            exactP = .false.
            print 500, L, Pobs(L), Pstd(L), x
          endif
          Pobs(L) = Pstd(L)    ! make it as accurate as possible
        enddo
         
        if( .NOT.exactP ) then
          print 510
          i=8
          call softexit('INIT_RTA',i)
        endif

C       check # coefs correct for sets 1-7 ---\
        if( N1COEF.ne.28 ) then              ! 4+8+11+5
          print 610, 'N1COEF', N1COEF
          i=8
          call softexit('INIT_RTA',i)
        endif

        if( N2COEF.ne.33 ) then     ! 4+8+10+11
          print 610, 'N2COEF', N2COEF 
          i=8
          call softexit('INIT_RTA',i)
        endif

        if( N3COEF.ne.32 ) then     ! 4+8+9+11
          print 610, 'N3COEF', N3COEF
          i=8
          call softexit('INIT_RTA',i)
        endif

        if( N4COEF.ne.42 ) then     ! 4+11+11+3+13
          print 610, 'N4COEF', N4COEF
          i=8
          call softexit('INIT_RTA',i)
        endif

        if( N5COEF.ne.19 ) then     ! 4+11+3+1
          print 610, 'N5COEF', N5COEF
          i=8
          call softexit('INIT_RTA',i)
        endif
       
        if( N6COEF.ne.20 ) then    ! 4+8+7+1
          print 610, 'N6COEF', N6COEF
          i=8
          call softexit('INIT_RTA',i)
        endif
             
        if( N7COEF.ne.26 ) then  ! 4+8+13+1
          print 610, 'N7COEF', N7COEF 
          i=8
          call softexit('INIT_RTA',i)
        endif

C       check # coefs correct for co2
        if( NCO2.ne.4 ) then
          print 610, 'NCO2', NCO2
          i=8
          call softexit('INIT_RTA',i)
        endif

C       OPTRAN water
        if( NOWAVG.ne.4 ) then       
          print 610, 'NOWAVG', NOWAVG  
          i=8  
          call softexit('INIT_RTA',i)
        endif

        if( NH2O.ne.9 ) then
          print 610, 'NH2O', NH2O
          i=8
          call softexit('INIT_RTA',i)
        endif

      endif

      firstcall = .false.

C     *********************************************************************
C                    EXECUTABLE CODE
C     *********************************************************************
C
C     Initialize the sum terms to zero

      PNORM=0.0E+0
      TZ=0.0E+0
      WZREF=0.0E+0
      WZ=0.0E+0
      XZREF=0.0E+0
      XZ=0.0E+0
      OZREF=0.0E+0
      OZ=0.0E+0
      TOZ=0.0E+0
      CZREF=0.0E+0
      CZ=0.0E+0
      MZREF=0.0E+0
      MZ=0.0E+0
      TMZ=0.0E+0

C     --------------------
C     Loop over the layers
C     --------------------

c     NOTE:  for NAST (LTOP>1),the  predictors are set to zero in init_com()

      do 8000 L = LTOP,NLAY
C
C       ---------------------------
C       Calculate the basic profile
C       dependent predictors.
C         ---------------------------
C
        IF (L .EQ. LTOP) THEN
          dz = 1.0E+5*(ALT(L) - ALT(L+1))  ! close enough
          PDP=Patm(L)*( Patm(L+1) - Patm(L))   ! NOTE: this works in mB or atmospheres
          TRZ=0.0E+0
          TAZ_O=0.0E+0
          TAZ_M=0.0E+0
        ELSE
          dz = 1.0E+5*(ALT(L-1) - ALT(L))  ! thickness of layer in cm's
          PDP=Patm(L)*( Patm(L) - Patm(L-1) )
          PNORM=PNORM + PDP

C         Note: TRZ, TOZ, and TMZ use layer-above terms
          TZ=TZ + PDP*TR
          TRZ=TZ/PNORM

          TOZ=TOZ + PDP*DT*A_O
          TAZ_O=TOZ/PNORM

          TMZ=TMZ + PDP*TR*A_M
          TAZ_M=TMZ/PNORM

        ENDIF
C
C       Temperature terms
        DT = Tavg(L) - Tstd(L)
        TR = Tavg(L)/Tstd(L)
C
C       Fixed term
c       A_F=PFAMNT(L)/RFAMNT(L)
        Pwater = (wcd(L)*Tavg(L))/(dz*273.15*2.6867E+19)
        A_F = (1.0 - Pmult*Pwater/Patm(L))/FX(L)
C
C       Water terms
        A_W=wcd(L)/wcdstd(L)         ! W
        WZREF=WZREF + PDP*wcdstd(L) 
        WZ=WZ + PDP*wcd(L)
        AZ_W=WZ/WZREF                ! WZ
C
C       Ozone terms
        A_O=ocd(L)/ocdstd(L)         ! O
        XZREF=XZREF + ocdstd(L)
        XZ=XZ + ocd(L)
        XZ_O=XZ/XZREF                ! Ox
        OZREF=OZREF + PDP*ocdstd(L)
        OZ=OZ + PDP*ocd(L)
        AZ_O=OZ/OZREF                ! Oz

        if(init_traceflg) then
          CPRED(9,L) = cocd(L)/costd(L)   ! A_C(L)  ... save for later
          CZREF=CZREF+PDP*costd(L)
          CZ=CZ + PDP*cocd(L)    !!!!! cocd
          CPRED(12,L) = CZ/CZREF          ! AZ_C(L) ... save for later

          MPRED(10,L) = ch4cd(L)/ch4std(L) ! A_M(L)  ... save for later
          MZREF=MZREF + PDP*ch4std(L)
          MZ=MZ + PDP*ch4cd(L)
          MPRED(11,L) = MZ/MZREF           ! AZ_M(L) ... save for later
        endif

C       --------------------------------------------
C       Load up the local variables for a = SECANG()
C       --------------------------------------------

        a = SECANG(L)                  ! a
        a2 = a*a                       ! a^^2

        TJUNKS=TR*TR                   ! Tr^2
        FIXMUL(L)=A_F

        OJUNKA=a*A_O                   ! O*a
        OJUNKR=SQRT( OJUNKA )          ! sqrt(O*a)
        OJUNKZ=OJUNKA/XZ_O             ! O*a/Ox
        OJUNKX=a*XZ_O                  ! a*Ox
 
        WJUNKA=a*A_W                   ! W*a
        WJUNKR=SQRT( WJUNKA )          ! sqrt(W*a)
        WJUNKS=WJUNKA*WJUNKA           ! (W*a)^2
        WJUNKZ=WJUNKA*A_W/AZ_W         ! W*a*W/Wz
        WJUNK4=SQRT( WJUNKR )          ! root^4(W*a)

        A_M  = MPRED(10,L)             ! M
        AZ_M = MPRED(11,L)             ! Mz
        MJUNKA=a*A_M                   ! M*a
        MJUNKR=SQRT(MJUNKA)            ! sqrt(M*a)
        MJUNKZ=a*AZ_M                  ! a*Mz

        A_C  = CPRED( 9,L)             ! C
        AZ_C = CPRED(12,L)             ! Cz
        CJUNKA=a*A_C                   ! C*a
        CJUNKR=SQRT( CJUNKA )          ! sqrt(C*a)
        CJUNKS=CJUNKA*CJUNKA           ! (C*a)^2
        CJUNKZ=CJUNKA*A_C/AZ_C         ! (C*a)*C/Cz

C       ---------------
C       Water continuum (for FWO, FOW, FMW, FCOW)
C       ---------------
        CONPRD(1,L)=WJUNKA/TJUNKS          ! a*W/Tr^2
        CONPRD(2,L)=CONPRD(1,L)*A_W/TJUNKS ! a*(W/Tr^2)^2
        CONPRD(3,L)=WJUNKA/TR              ! a*W/Tr
        CONPRD(4,L)=CONPRD(3,L)*A_W        ! a*W^2/Tr

        FPRED( 1,L) = a                    ! a
        FPRED( 2,L) = a2                   ! a*a
        FPRED( 3,L) = a*TR                 ! a*Tr
        FPRED( 4,L) = a*TJUNKS             ! a*Tr^2
        FPRED( 5,L) = TR                   ! Tr
        FPRED( 6,L) = TJUNKS               ! Tr**2
        FPRED( 7,L) = a*TRZ                ! a*Trz
        FPRED( 8,L) = a*TRZ/TR             ! a*Trz/Tr 
        FPRED( 9,L) = a2*TR                ! a^2*Tr
        FPRED(10,L) = a2*a                 ! a^3
        FPRED(11,L) = SQRT(a)              ! sqrt(a)
        FPRED(12,L) = a2*TRZ               ! a^2*Trz
        FPRED(13,L) = TRZ                  ! Trz
        FPRED(14,L) = sqrt(TR)             ! sqrt(Tr)

C       -----
C       Ozone
C       -----
C
C       Ozone predictors for FWO = set1
        OPRED( 1,L)=OJUNKA                 ! O*a 
        OPRED( 2,L)=OJUNKR                 ! sqrt(O*a) 
        OPRED( 3,L)=OJUNKA*DT              ! O*a*dT
        OPRED( 4,L)=OJUNKA*OJUNKA          ! (O*a)^2
        OPRED( 5,L)=OJUNKR*DT              ! sqrt(O*a)*dT
        OPRED( 6,L)=OJUNKZ*A_O             ! O*a*O/Ox
        OPRED( 7,L)=OJUNKR*A_O/XZ_O        ! sqrt(O*a)*O/Ox
        OPRED( 8,L)=OJUNKZ*AZ_O            ! O*a*Oz/Ox
        OPRED( 9,L)=OJUNKA*SQRT( OJUNKX )  ! O*a*sqrt(Ox*a)
        OPRED(10,L)=OJUNKA*TAZ_O*a         ! O*a*TOz*a
        OPRED(11,L)=OJUNKA*OJUNKR          ! (O*a)^3/2
        OPRED(12,L)=OJUNKA*a               ! O*a^2
        OPRED(13,L)=OJUNKA*OJUNKZ*A_O      ! (O*a)^2*O/Ox
        OPRED(14,L)=OJUNKA**3              ! (O*a)^3
        OPRED(15,L)=A_O                    ! O

C       -----
C       Water
C       -----
C
C       Water predictors for all sets
        WPRED( 1,L)=WJUNKA                ! W*a
        WPRED( 2,L)=WJUNKR                ! sqrt(W*a)
        WPRED( 3,L)=WJUNKZ                ! W*a*W/Wz 
        WPRED( 4,L)=WJUNKA*DT             ! W*a*dT
        WPRED( 5,L)=WJUNKS                ! (W*a)^2
        WPRED( 6,L)=WJUNKR*DT             ! sqrt(W*a)*dT
        WPRED( 7,L)=WJUNK4                ! root^4(W*a)
        WPRED( 8,L)=WJUNKZ/WJUNKR         ! sqrt(W*a)*W/Wz
        WPRED( 9,L)=WJUNKS*WJUNKA         ! (W*a)^3
        WPRED(10,L)=A_W                   ! W
        WPRED(11,L)=WJUNKA*DT*ABS( DT )   ! W*a*dT*|dT|
        WPRED(12,L)=WJUNKA*OJUNKX         ! W*a*Ox*a
        WPRED(13,L)=WJUNKA*OJUNKX*OJUNKX  ! W*a*(Ox*a)^2
        WPRED(14,L)=WJUNKR*MJUNKZ         ! W*a*Mz*a 
        WPRED(15,L)=WJUNKA*a              ! W*a^2 
        WPRED(16,L)=WJUNKA*AZ_C*a         ! W*a*Cz*a
        WPRED(17,L)=WJUNKA*DT*a           ! W*a^2*dT
        WPRED(18,L)=WJUNKA*WJUNKR         ! (W*a)^3/2
        WPRED(19,L)=WJUNKA*WJUNKR*DT      ! (W*a)^3/2*dT
        WPRED(20,L)=WJUNKZ*WJUNKR         ! (W*a)^3/2*W/Wz = (W*a*W/Wz)*sqrt(W*a) 
        WPRED(21,L)=WJUNKA*WJUNK4         ! (W*a)^5/4
        WPRED(22,L)=WJUNKA*WJUNKZ         ! (W*a)^2*W/Wz
        WPRED(23,L)=WJUNKA*A_W            ! W^2*a
        WPRED(24,L)=WJUNKS/WJUNK4         ! (W*a)^7/4
        WPRED(25,L)=WJUNKA/WJUNK4         ! (W*a)^3/4

        if(init_traceflg) then
C         -------
C         Methane for FMW = set3
C         -------

          MPRED(1,L)=MJUNKA               ! M*a
          MPRED(2,L)=MJUNKR               ! sqrt(M*a)
C--       MPRED(3,L)=MJUNKA*DT            ! M*a*dT   == MPRED(1,L)*DT 
          MPRED(4,L)=MJUNKA*MJUNKA        ! (M*a)^2
          MPRED(5,L)=MJUNKA*a             ! M*a^2
          MPRED(6,L)=MJUNKZ               ! Mz*a
C--       MPRED(7,L)=A_M*DT               ! M*dT     == MPRED(10,L)*DT
          MPRED(8,L)=TAZ_M*a              ! TMz*a 
          MPRED(9,L)=SQRT( MJUNKZ )       ! sqrt(Mz*a)
C--       MPRED(10,L) = A_M               ! M   saved above, used in MPRED(7,L) computation
C--       MPRED(11,L) = AZ_M              ! Mz  saved above

C         ---------------
C         Carbon monoxide for FCOW = set4
C         ---------------

          CPRED(1,L)=CJUNKA              ! C*a 
          CPRED(2,L)=CJUNKR              ! sqrt(C*a)
C--       CPRED(3,L)=CJUNKA*DT           ! C*a*dT     == CPRED(1,L)*DT  --- below
          CPRED(4,L)=CJUNKS              ! (C*a)^2 
          CPRED(5,L)=CJUNKZ              ! C*a*C/Cz
C--       CPRED(6,L)=CJUNKR*DT           ! sqrt(C*a)*dT == CPRED(2,L)*DT  --- below
          CPRED(7,L)=SQRT( CJUNKR )      ! root^4(C*a)
          CPRED(8,L)=CJUNKZ/CJUNKR       ! sqrt(C*a)*C/Cz
C--       CPRED(9,L) = A_C               ! C   done above to save A_C(L)
          CPRED(10,L)=CJUNKA*a           ! a^2*C
          CPRED(11,L)=CJUNKR*a           ! C*sqrt(C*a)
C--       CPRED(12,L) = AZ_C             ! Cz  done above to save AZ_C(L)

C         ---------------
C         CO2 perturbation coefs
C         ---------------
            CO2PRD(1,L)=a                ! a

        endif

C       temperature dependent methane terms
C       -----------------------------------
        MPRED(3,L) = MPRED(1,L)*DT       ! M*a*dT  ==  MJUNKA*DT
        MPRED(7,L) = MPRED(10,L)*DT      ! M*dT    == A_M*DT

C       temperature dependent carbon monoxide terms
C       -------------------------------------------
 
        CPRED(3,L)=CPRED(1,L)*DT         ! C*a*dT
        CPRED(6,L)=CPRED(2,L)*DT         ! sqrt(C*a)*dT

C       temperature dependent carbon dioxide terms
C       ------------------------------------------
        CO2PRD(2,L)=TR                   ! Tr
        CO2PRD(3,L)=a*TR                 ! a*Tr
        CO2PRD(4,L)=a*TJUNKS             ! a*Tr^2

C       ####################################################
C       ########## NOW initialize the SOLAR TERMS ##########
C       ####################################################

        if(solzang.ge.solzang_limit) goto 8000   ! cossun() = 0.0, TAUZSN=0.0 (ir_rta)

        a = SECSUN(L)                  ! a
        a2 = a*a                       ! a^2

        OJUNKA=a*A_O                   ! O*a
        OJUNKR=SQRT( OJUNKA )          ! sqrt(O*a)
        OJUNKZ=OJUNKA/XZ_O             ! O*a/Ox
        OJUNKX=a*XZ_O                  ! a*Ox

        WJUNKA=a*A_W                   ! W*a
        WJUNKR=SQRT( WJUNKA )          ! sqrt(W*a)
        WJUNKS=WJUNKA*WJUNKA           ! (W*a)^2
        WJUNKZ=WJUNKA*A_W/AZ_W         ! W*a*W/Wz
        WJUNK4=SQRT( WJUNKR )          ! root^4(W*a)

        MJUNKZ=a*AZ_M                  ! a*Mz  , also used in water

        CJUNKA=a*A_C                   ! C*a
        CJUNKR=SQRT( CJUNKA )          ! SQRT(C*a)
        CJUNKS=CJUNKA*CJUNKA           ! (C*a)^2
        CJUNKZ=CJUNKA*A_C/AZ_C         ! (C*a)*C/Cz

        CONPRD_SOL(1,L)=WJUNKA/TJUNKS               ! a*W/Tr^2
        CONPRD_SOL(2,L)=CONPRD_SOL(1,L)*A_W/TJUNKS  ! a*(W/Tr^2)^2
        CONPRD_SOL(3,L)=WJUNKA/TR                   ! a*W/Tr
        CONPRD_SOL(4,L)=CONPRD_SOL(3,L)*A_W         ! a*W^2/Tr

        FPRED_SOL( 1,L) = a            ! a
        FPRED_SOL( 2,L) = a2           ! a^2
        FPRED_SOL( 3,L) = a*TR         ! a*Tr
        FPRED_SOL( 4,L) = a*TJUNKS     ! a*Tr^2
        FPRED_SOL( 5,L) = TR           ! Tr
        FPRED_SOL( 6,L) = TJUNKS       ! Tr**2
        FPRED_SOL( 7,L) = a*TRZ        ! a*Trz
        FPRED_SOL( 8,L) = a*TRZ/TR     ! a*Trz/Tr 
        FPRED_SOL( 9,L) = a2*TR        ! a^2*Tr
        FPRED_SOL(10,L) = a2*a         ! a^3
        FPRED_SOL(11,L) = SQRT(a)      ! sqrt(a)
        FPRED_SOL(12,L) = a2*TRZ       ! a^2*Trz
        FPRED_SOL(13,L) = TRZ          ! Trz
        FPRED_SOL(14,L) = sqrt(TR)     ! sqrt(Tr)

C       -----
C       Ozone
C       -----

        OPRED_SOL( 1,L)=OJUNKA               ! O*a 
        OPRED_SOL( 2,L)=OJUNKR               ! sqrt(O*a) 
        OPRED_SOL( 3,L)=OJUNKA*DT            ! O*a*dT
        OPRED_SOL( 4,L)=OJUNKA*OJUNKA        ! (O*a)^2
        OPRED_SOL( 5,L)=OJUNKR*DT            ! sqrt(O*a)*dT
        OPRED_SOL( 6,L)=OJUNKZ*A_O           ! O*a*O/Ox
        OPRED_SOL( 7,L)=OJUNKR*A_O/XZ_O      ! sqrt(O*a)*O/Ox
        OPRED_SOL( 8,L)=OJUNKZ*AZ_O          ! O*a*Oz/Ox
        OPRED_SOL( 9,L)=OJUNKA*SQRT(OJUNKX)  ! O*a*sqrt(Ox*a)
        OPRED_SOL(10,L)=OJUNKA*TAZ_O*a       ! O*a*TOz*a
        OPRED_SOL(11,L)=OJUNKA*OJUNKR        ! (O*a)^3/2
        OPRED_SOL(12,L)=OJUNKA*a             ! O*a^2
        OPRED_SOL(13,L)=OJUNKA*OJUNKZ*A_O    ! (O*a)^2*O/Ox
        OPRED_SOL(14,L)=OJUNKA**3            ! (O*a)^3
        OPRED_SOL(15,L)=A_O                  ! O

C       -----
C       Water
C       -----

        WPRED_SOL( 1,L)=WJUNKA               ! W*a
        WPRED_SOL( 2,L)=WJUNKR               ! sqrt(W*a)
        WPRED_SOL( 3,L)=WJUNKZ               ! W*a*W/Wz 
        WPRED_SOL( 4,L)=WJUNKA*DT            ! W*a*dT
        WPRED_SOL( 5,L)=WJUNKS               ! (W*a)^2
        WPRED_SOL( 6,L)=WJUNKR*DT            ! sqrt(W*a)*dT
        WPRED_SOL( 7,L)=WJUNK4               ! root^4(W*a)
        WPRED_SOL( 8,L)=WJUNKZ/WJUNKR        ! sqrt(W*a)*W/Wz
        WPRED_SOL( 9,L)=WJUNKS*WJUNKA        ! (W*a)^3
        WPRED_SOL(10,L)=A_W                  ! W
        WPRED_SOL(11,L)=WJUNKA*DT*ABS(DT)    ! W*a*dT*|dT|
        WPRED_SOL(12,L)=WJUNKA*OJUNKX        ! W*a*Ox*a
        WPRED_SOL(13,L)=WJUNKA*OJUNKX*OJUNKX ! W*a*(Ox*a)^2
        WPRED_SOL(14,L)=WJUNKR*MJUNKZ        ! W*a*Mz*a 
        WPRED_SOL(15,L)=WJUNKA*a             ! W*a^2 
        WPRED_SOL(16,L)=WJUNKA*AZ_C*a        ! W*a*Cz*a
        WPRED_SOL(17,L)=WJUNKA*DT*a          ! W*a^2*dT
        WPRED_SOL(18,L)=WJUNKA*WJUNKR        ! (W*a)^3/2
        WPRED_SOL(19,L)=WJUNKA*WJUNKR*DT     ! (W*a)^3/2*dT
        WPRED_SOL(20,L)=WJUNKZ*WJUNKR        ! (W*a)^3/2*W/Wz
        WPRED_SOL(21,L)=WJUNKA*WJUNK4        ! (W*a)^5/4
        WPRED_SOL(22,L)=WJUNKA*WJUNKZ        ! (W*a)^2*W/Wz
        WPRED_SOL(23,L)=WJUNKA*A_W           ! W^2*a
        WPRED_SOL(24,L)=WJUNKS/WJUNK4        ! (W*a)^7/4
        WPRED_SOL(25,L)=WJUNKA/WJUNK4        ! (W*a)^3/4

        if(init_traceflg) then

C         ---------------
C         Carbon monoxide for FCOW = set4 SOLAR TERMS
C         ---------------
C        
          CPRED_SOL(1,L)=CJUNKA
          CPRED_SOL(2,L)=CJUNKR
C--       CPRED_SOL(3,L)=CJUNKA*DT == CPRED_SOL(1,L)*DT  SOLAR
          CPRED_SOL(4,L)=CJUNKS   
          CPRED_SOL(5,L)=CJUNKZ
C--       CPRED_SOL(6,L)=CJUNKR*DT == CPRED_SOL(2,L)*DT  SOLAR
          CPRED_SOL(7,L)=SQRT( CJUNKR )
          CPRED_SOL(8,L)=CJUNKZ/CJUNKR
C--       CPRED_SOL(9,L)=A_C         ! done above to save AZ_C
          CPRED_SOL(10,L)=CJUNKA*a
          CPRED_SOL(11,L)=CJUNKR*a
C--       CPRED_SOL(12,L)=AZ_C       ! done above to save AZ_C


C         ---------------
C         CO2 perturbation coefs
C         ---------------
          CO2PRD_SOL(1,L)=a

        endif

        CPRED_SOL(3,L)=CPRED_SOL(1,L)*DT
        CPRED_SOL(6,L)=CPRED_SOL(2,L)*DT

        CO2PRD_SOL(2,L)=TR
        CO2PRD_SOL(3,L)=a*TR
        CO2PRD_SOL(4,L)=a*TJUNKS

 8000 continue
      init_traceflg = .false.   ! done for now

C     ####################################################
C     ######### NOW initialize the OPTRAN TERMS ##########
C     ####################################################

      if(.NOT.OPTRAN) goto 9999


C-- INIT_OPTRAN/CALOKW -------\

C- For AIRS ONLY!...taken from USEFAST
C-- From UMBC-usefast -> should it go here?  ---\
C-- Calc OPTRAN absorption coefficient scaling factor ---\
       
C      Initialize amount above sums
       WAZSUM=0.0E+0
       WPZSUM=0.0E+0
       WTZSUM=0.0E+0
C
C      ---------------------------------------
C      Calculate raw predictors for all layers
C      ---------------------------------------
       DO L = LTOP,NLAY
C
C         Layer amount*angle
          WAANG(L)=SECANG(L)* wcd(L)*kmole
C
C         Center-of-layer-to-space amount*angle
C         Note: do this before updating AZSUM
          WAZ(L)=5.0E-1*WAANG(L) + WAZSUM
C
C         Bottom-of-layer-to-space amount sum
          WAZSUM = WAZSUM + WAANG(L)
C
C         Pressure above sum
c         WPZSUM=WAANG(L)*P(L) + WPZSUM
          WPZSUM=WAANG(L)*Patm(L) + WPZSUM
          PZ(L)=WPZSUM/WAZSUM
C
C         Temperature above sum
c         WTZSUM=WAANG(L)*T(L) + WTZSUM
          WTZSUM=WAANG(L)*Tavg(L) + WTZSUM
          TZopt(L)=WTZSUM/WAZSUM
C         
       ENDDO
C
C      --------------------------------------------------
C      Find the max OPTRAN level that is less than WAZ(1)
C      --------------------------------------------------
       LOPMIN=1
 30    IF (WAZOP(LOPMIN+1) .LT. WAZ(LTOP)) THEN
          LOPMIN=LOPMIN + 1
          GOTO 30
       ENDIF
C
C      Initialize the upper and lower (pressure) layer index
       LL=LTOP
       LU=LTOP+1
       LAST=.FALSE.
C         
C      ----------------------------------------
C      Loop over the OPTRAN layers (while loop)
C      ----------------------------------------
       LOP=LOPMIN
 10    IF (LOP .LE. MXOWLY) THEN
C
C         --------------------------------------------------------
C         Find the two pressure layers closest to the OPTRAN layer
C         --------------------------------------------------------
 20       IF (WAZ(LU) .LT. WAZOP(LOP)) THEN
             IF (LU .LT. NLAY) THEN
                LL=LU
                LU=LU + 1
                GOTO 20
             ELSE  
                LAST=.TRUE.
             ENDIF
          ENDIF

          if(LU.gt.NLAY) then
            print 21, LU, NLAY
            i = 8
            call softexit('INIT_RTA',i)
          endif

C         Compute the interpolation factor
          DA=(WAZOP(LOP) - WAZ(LL))/(WAZ(LU) - WAZ(LL))
C
C         Do the interpolation
c          POP=( DA*(  P(LU) -  P(LL) ) +  P(LL) )/WAVGOP(1,LOP)  
c          TOP=( DA*(  T(LU) -  T(LL) ) +  T(LL) )/WAVGOP(2,LOP)
           POP=( DA*( Patm(LU)-Patm(LL) ) + Patm(LL) )/WAVGOP(1,LOP)
           TOP=( DA*(Tavg(LU)-Tavg(LL) ) +  Tavg(LL) )/WAVGOP(2,LOP)
          PZOP=( DA*( PZ(LU) - PZ(LL) ) + PZ(LL) )/WAVGOP(3,LOP)
          TZOP=( DA*(TZopt(LU)-TZopt(LL))+TZopt(LL) )/WAVGOP(4,LOP)
          ANGOP=DA*( SECANG(LU) - SECANG(LL) ) + SECANG(LL)
C
C         Assign the predictors
          H2OPRD(1,LOP)=1.0E+0
          H2OPRD(2,LOP)=POP
          H2OPRD(3,LOP)=TOP
          H2OPRD(4,LOP)=SQRT( POP )
          H2OPRD(5,LOP)=TOP**2
          H2OPRD(6,LOP)=POP*TOP
          H2OPRD(7,LOP)=ANGOP 
          H2OPRD(8,LOP)=PZOP
          H2OPRD(9,LOP)=TZOP
C
C         Update LOP and loop
          IF (LAST) THEN
             LOPMAX=LOP
C            Set LOP > MXOWLY to exit loop over LOP
             LOP=MXOWLY + 1
          ELSE
             LOP=LOP + 1
          ENDIF
          GOTO 10
C
       ENDIF
C      End while loop over LOP 
C
C      -----------------
C      Initialize LOPUSE
C      -----------------
       DO LOP=1,MXOWLY
          LOPUSE(LOP)=.FALSE.
       ENDDO
C            
C      ---------------------------------------
C      Determine what OPTRAN layers are needed
C      ---------------------------------------
C      Initialize LOPL and LOPU
       LOPL=LOPMIN
       LOPU=LOPMIN + 1
C
C      Loop over the AIRS pressure layers
       DO L = LTOP, NLAY
C         Find the two OPTRAN levels that bracket the AIRS layer
 40       IF (WAZOP(LOPU) .LT. WAZ(L) .AND. LOPU .LT. LOPMAX) THEN
             LOPL=LOPU  
             LOPU=LOPU + 1
             GOTO 40
          ENDIF
C            
          LOPUSE(LOPL)=.TRUE.
          LOPUSE(LOPU)=.TRUE.
C         Assign the lower OPTRAN level
          LOPLOW(L)=LOPL
C         Assign the interpolation fraction
          DAOP(L)=(WAZ(L) - WAZOP(LOPL))/(WAZOP(LOPU) - WAZOP(LOPL))
      ENDDO

 9999 RETURN
  21  format('init_rta: OPTRAN layer=',i4,' >',i4)
 100  FORMAT('init_rta: physics not supported, iphys=',i3)
 400  FORMAT('init_rta: SEVERE ERROR NLAY=',i3,' > ',
     +  i3,'=ntrlev_ir',' Psurf=',f9.2,' bottom of RTA =',f9.2)
 500  FORMAT(5x,'L=',i3,' Pobs(L)=',f10.4,' Pstd(L)=',f10.4,
     +   ' diff=',f10.7)
 510  FORMAT('init_rta: fatal error, pressure scales not equal')
 610  FORMAT('init_rta: ',a6,'=',i3,' which is wrong.')
      END
