LCOV - code coverage report
Current view: top level - physics/rrtmg/aer_src - rrtmg_lw_taumol.f90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 1771 1855 95.5 %
Date: 2025-03-13 18:55:17 Functions: 17 17 100.0 %

          Line data    Source code
       1             : !     path:      $Source: /storm/rc1/cvsroot/rc/rrtmg_lw/src/rrtmg_lw_taumol.f90,v $
       2             : !     author:    $Author: mike $
       3             : !     revision:  $Revision: 1.7 $
       4             : !     created:   $Date: 2009/10/20 15:08:37 $
       5             : !
       6             :       module rrtmg_lw_taumol
       7             : 
       8             : !  --------------------------------------------------------------------------
       9             : ! |                                                                          |
      10             : ! |  Copyright 2002-2009, Atmospheric & Environmental Research, Inc. (AER).  |
      11             : ! |  This software may be used, copied, or redistributed as long as it is    |
      12             : ! |  not sold and this copyright notice is reproduced on each copy made.     |
      13             : ! |  This model is provided as is without any express or implied warranties. |
      14             : ! |                       (http://www.rtweb.aer.com/)                        |
      15             : ! |                                                                          |
      16             : !  --------------------------------------------------------------------------
      17             : 
      18             : ! ------- Modules -------
      19             : 
      20             :       use shr_kind_mod, only: r8 => shr_kind_r8
      21             : 
      22             :       use parrrtm,  only: nbndlw, maxxsec, ngptlw
      23             :       use rrlw_con, only: oneminus
      24             :       use rrlw_wvn, only: nspa, nspb
      25             : 
      26             :       implicit none
      27             : 
      28             :       contains
      29             : 
      30             : !----------------------------------------------------------------------------
      31      486000 :       subroutine taumol(nlayers, pavel, wx, coldry, &
      32      486000 :                         laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
      33      486000 :                         colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
      34      486000 :                         colbrd, fac00, fac01, fac10, fac11, &
      35      486000 :                         rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
      36      486000 :                         rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
      37      486000 :                         rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
      38      486000 :                         selffac, selffrac, indself, forfac, forfrac, indfor, &
      39      486000 :                         minorfrac, scaleminor, scaleminorn2, indminor, &
      40      486000 :                         fracs, taug)
      41             : !----------------------------------------------------------------------------
      42             : 
      43             : ! *******************************************************************************
      44             : ! *                                                                             *
      45             : ! *                  Optical depths developed for the                           *
      46             : ! *                                                                             *
      47             : ! *                RAPID RADIATIVE TRANSFER MODEL (RRTM)                        *
      48             : ! *                                                                             *
      49             : ! *                                                                             *
      50             : ! *            ATMOSPHERIC AND ENVIRONMENTAL RESEARCH, INC.                     *
      51             : ! *                        131 HARTWELL AVENUE                                  *
      52             : ! *                        LEXINGTON, MA 02421                                  *
      53             : ! *                                                                             *
      54             : ! *                                                                             *
      55             : ! *                           ELI J. MLAWER                                     * 
      56             : ! *                         JENNIFER DELAMERE                                   * 
      57             : ! *                         STEVEN J. TAUBMAN                                   *
      58             : ! *                         SHEPARD A. CLOUGH                                   *
      59             : ! *                                                                             *
      60             : ! *                                                                             *
      61             : ! *                                                                             *
      62             : ! *                                                                             *
      63             : ! *                       email:  mlawer@aer.com                                *
      64             : ! *                       email:  jdelamer@aer.com                              *
      65             : ! *                                                                             *
      66             : ! *        The authors wish to acknowledge the contributions of the             *
      67             : ! *        following people:  Karen Cady-Pereira, Patrick D. Brown,             *  
      68             : ! *        Michael J. Iacono, Ronald E. Farren, Luke Chen, Robert Bergstrom.    *
      69             : ! *                                                                             *
      70             : ! *******************************************************************************
      71             : ! *                                                                             *
      72             : ! *  Revision for g-point reduction: Michael J. Iacono, AER, Inc.               *
      73             : ! *                                                                             *
      74             : ! *******************************************************************************
      75             : ! *     TAUMOL                                                                  *
      76             : ! *                                                                             *
      77             : ! *     This file contains the subroutines TAUGBn (where n goes from            *
      78             : ! *     1 to 16).  TAUGBn calculates the optical depths and Planck fractions    *
      79             : ! *     per g-value and layer for band n.                                       *
      80             : ! *                                                                             *
      81             : ! *  Output:  optical depths (unitless)                                         *
      82             : ! *           fractions needed to compute Planck functions at every layer       *
      83             : ! *               and g-value                                                   *
      84             : ! *                                                                             *
      85             : ! *     COMMON /TAUGCOM/  TAUG(MXLAY,MG)                                        *
      86             : ! *     COMMON /PLANKG/   FRACS(MXLAY,MG)                                       *
      87             : ! *                                                                             *
      88             : ! *  Input                                                                      *
      89             : ! *                                                                             *
      90             : ! *     COMMON /FEATURES/ NG(NBANDS),NSPA(NBANDS),NSPB(NBANDS)                  *
      91             : ! *     COMMON /PRECISE/  ONEMINUS                                              *
      92             : ! *     COMMON /PROFILE/  NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),                    *
      93             : ! *     &                 PZ(0:MXLAY),TZ(0:MXLAY)                               *
      94             : ! *     COMMON /PROFDATA/ LAYTROP,                                              *
      95             : ! *    &                  COLH2O(MXLAY),COLCO2(MXLAY),COLO3(MXLAY),             *
      96             : ! *    &                  COLN2O(MXLAY),COLCO(MXLAY),COLCH4(MXLAY),             *
      97             : ! *    &                  COLO2(MXLAY)
      98             : ! *     COMMON /INTFAC/   FAC00(MXLAY),FAC01(MXLAY),                            *
      99             : ! *    &                  FAC10(MXLAY),FAC11(MXLAY)                             *
     100             : ! *     COMMON /INTIND/   JP(MXLAY),JT(MXLAY),JT1(MXLAY)                        *
     101             : ! *     COMMON /SELF/     SELFFAC(MXLAY), SELFFRAC(MXLAY), INDSELF(MXLAY)       *
     102             : ! *                                                                             *
     103             : ! *     Description:                                                            *
     104             : ! *     NG(IBAND) - number of g-values in band IBAND                            *
     105             : ! *     NSPA(IBAND) - for the lower atmosphere, the number of reference         *
     106             : ! *                   atmospheres that are stored for band IBAND per            *
     107             : ! *                   pressure level and temperature.  Each of these            *
     108             : ! *                   atmospheres has different relative amounts of the         *
     109             : ! *                   key species for the band (i.e. different binary           *
     110             : ! *                   species parameters).                                      *
     111             : ! *     NSPB(IBAND) - same for upper atmosphere                                 *
     112             : ! *     ONEMINUS - since problems are caused in some cases by interpolation     *
     113             : ! *                parameters equal to or greater than 1, for these cases       *
     114             : ! *                these parameters are set to this value, slightly < 1.        *
     115             : ! *     PAVEL - layer pressures (mb)                                            *
     116             : ! *     TAVEL - layer temperatures (degrees K)                                  *
     117             : ! *     PZ - level pressures (mb)                                               *
     118             : ! *     TZ - level temperatures (degrees K)                                     *
     119             : ! *     LAYTROP - layer at which switch is made from one combination of         *
     120             : ! *               key species to another                                        *
     121             : ! *     COLH2O, COLCO2, COLO3, COLN2O, COLCH4 - column amounts of water         *
     122             : ! *               vapor,carbon dioxide, ozone, nitrous ozide, methane,          *
     123             : ! *               respectively (molecules/cm**2)                                *
     124             : ! *     FACij(LAY) - for layer LAY, these are factors that are needed to        *
     125             : ! *                  compute the interpolation factors that multiply the        *
     126             : ! *                  appropriate reference k-values.  A value of 0 (1) for      *
     127             : ! *                  i,j indicates that the corresponding factor multiplies     *
     128             : ! *                  reference k-value for the lower (higher) of the two        *
     129             : ! *                  appropriate temperatures, and altitudes, respectively.     *
     130             : ! *     JP - the index of the lower (in altitude) of the two appropriate        *
     131             : ! *          reference pressure levels needed for interpolation                 *
     132             : ! *     JT, JT1 - the indices of the lower of the two appropriate reference     *
     133             : ! *               temperatures needed for interpolation (for pressure           *
     134             : ! *               levels JP and JP+1, respectively)                             *
     135             : ! *     SELFFAC - scale factor needed for water vapor self-continuum, equals    *
     136             : ! *               (water vapor density)/(atmospheric density at 296K and        *
     137             : ! *               1013 mb)                                                      *
     138             : ! *     SELFFRAC - factor needed for temperature interpolation of reference     *
     139             : ! *                water vapor self-continuum data                              *
     140             : ! *     INDSELF - index of the lower of the two appropriate reference           *
     141             : ! *               temperatures needed for the self-continuum interpolation      *
     142             : ! *     FORFAC  - scale factor needed for water vapor foreign-continuum.        *
     143             : ! *     FORFRAC - factor needed for temperature interpolation of reference      *
     144             : ! *                water vapor foreign-continuum data                           *
     145             : ! *     INDFOR  - index of the lower of the two appropriate reference           *
     146             : ! *               temperatures needed for the foreign-continuum interpolation   *
     147             : ! *                                                                             *
     148             : ! *  Data input                                                                 *
     149             : ! *     COMMON /Kn/ KA(NSPA(n),5,13,MG), KB(NSPB(n),5,13:59,MG), SELFREF(10,MG),*
     150             : ! *                 FORREF(4,MG), KA_M'MGAS', KB_M'MGAS'                        *
     151             : ! *        (note:  n is the band number,'MGAS' is the species name of the minor *
     152             : ! *         gas)                                                                *
     153             : ! *                                                                             *
     154             : ! *     Description:                                                            *
     155             : ! *     KA - k-values for low reference atmospheres (key-species only)          *
     156             : ! *          (units: cm**2/molecule)                                            *
     157             : ! *     KB - k-values for high reference atmospheres (key-species only)         *
     158             : ! *          (units: cm**2/molecule)                                            *
     159             : ! *     KA_M'MGAS' - k-values for low reference atmosphere minor species        *
     160             : ! *          (units: cm**2/molecule)                                            *
     161             : ! *     KB_M'MGAS' - k-values for high reference atmosphere minor species       *
     162             : ! *          (units: cm**2/molecule)                                            *
     163             : ! *     SELFREF - k-values for water vapor self-continuum for reference         *
     164             : ! *               atmospheres (used below LAYTROP)                              *
     165             : ! *               (units: cm**2/molecule)                                       *
     166             : ! *     FORREF  - k-values for water vapor foreign-continuum for reference      *
     167             : ! *               atmospheres (used below/above LAYTROP)                        *
     168             : ! *               (units: cm**2/molecule)                                       *
     169             : ! *                                                                             *
     170             : ! *     DIMENSION ABSA(65*NSPA(n),MG), ABSB(235*NSPB(n),MG)                     *
     171             : ! *     EQUIVALENCE (KA,ABSA),(KB,ABSB)                                         *
     172             : ! *                                                                             *
     173             : !*******************************************************************************
     174             : 
     175             : ! ------- Declarations -------
     176             : 
     177             : ! ----- Input -----
     178             :       integer, intent(in) :: nlayers         ! total number of layers
     179             :       real(kind=r8), intent(in) :: pavel(nlayers)           ! layer pressures (mb) 
     180             :                                                       !    Dimensions: (nlayers)
     181             :       real(kind=r8), intent(in) :: wx(maxxsec,nlayers)            ! cross-section amounts (mol/cm2)
     182             :                                                       !    Dimensions: (maxxsec,nlayers)
     183             :       real(kind=r8), intent(in) :: coldry(nlayers)          ! column amount (dry air)
     184             :                                                       !    Dimensions: (nlayers)
     185             : 
     186             :       integer, intent(in) :: laytrop         ! tropopause layer index
     187             :       integer, intent(in) :: jp(nlayers)           ! 
     188             :                                                       !    Dimensions: (nlayers)
     189             :       integer, intent(in) :: jt(nlayers)           !
     190             :                                                       !    Dimensions: (nlayers)
     191             :       integer, intent(in) :: jt1(nlayers)          !
     192             :                                                       !    Dimensions: (nlayers)
     193             :       real(kind=r8), intent(in) :: planklay(nlayers,nbndlw)      ! 
     194             :                                                       !    Dimensions: (nlayers,nbndlw)
     195             :       real(kind=r8), intent(in) :: planklev(0:nlayers,nbndlw)     ! 
     196             :                                                       !    Dimensions: (nlayers,nbndlw)
     197             :       real(kind=r8), intent(in) :: plankbnd(nbndlw)        ! 
     198             :                                                       !    Dimensions: (nbndlw)
     199             : 
     200             :       real(kind=r8), intent(in) :: colh2o(nlayers)          ! column amount (h2o)
     201             :                                                       !    Dimensions: (nlayers)
     202             :       real(kind=r8), intent(in) :: colco2(nlayers)          ! column amount (co2)
     203             :                                                       !    Dimensions: (nlayers)
     204             :       real(kind=r8), intent(in) :: colo3(nlayers)           ! column amount (o3)
     205             :                                                       !    Dimensions: (nlayers)
     206             :       real(kind=r8), intent(in) :: coln2o(nlayers)          ! column amount (n2o)
     207             :                                                       !    Dimensions: (nlayers)
     208             :       real(kind=r8), intent(in) :: colco(nlayers)           ! column amount (co)
     209             :                                                       !    Dimensions: (nlayers)
     210             :       real(kind=r8), intent(in) :: colch4(nlayers)          ! column amount (ch4)
     211             :                                                       !    Dimensions: (nlayers)
     212             :       real(kind=r8), intent(in) :: colo2(nlayers)           ! column amount (o2)
     213             :                                                       !    Dimensions: (nlayers)
     214             :       real(kind=r8), intent(in) :: colbrd(nlayers)          ! column amount (broadening gases)
     215             :                                                       !    Dimensions: (nlayers)
     216             : 
     217             :       integer, intent(in) :: indself(nlayers)
     218             :                                                       !    Dimensions: (nlayers)
     219             :       integer, intent(in) :: indfor(nlayers)
     220             :                                                       !    Dimensions: (nlayers)
     221             :       real(kind=r8), intent(in) :: selffac(nlayers)
     222             :                                                       !    Dimensions: (nlayers)
     223             :       real(kind=r8), intent(in) :: selffrac(nlayers)
     224             :                                                       !    Dimensions: (nlayers)
     225             :       real(kind=r8), intent(in) :: forfac(nlayers)
     226             :                                                       !    Dimensions: (nlayers)
     227             :       real(kind=r8), intent(in) :: forfrac(nlayers)
     228             :                                                       !    Dimensions: (nlayers)
     229             : 
     230             :       integer, intent(in) :: indminor(nlayers)
     231             :                                                       !    Dimensions: (nlayers)
     232             :       real(kind=r8), intent(in) :: minorfrac(nlayers)
     233             :                                                       !    Dimensions: (nlayers)
     234             :       real(kind=r8), intent(in) :: scaleminor(nlayers)
     235             :                                                       !    Dimensions: (nlayers)
     236             :       real(kind=r8), intent(in) :: scaleminorn2(nlayers)
     237             :                                                       !    Dimensions: (nlayers)
     238             : 
     239             :       real(kind=r8), dimension(nlayers), intent(in) :: &                  !
     240             :                        fac00, fac01, &          !    Dimensions: (nlayers)
     241             :                        fac10, fac11 
     242             :       real(kind=r8), dimension(nlayers), intent(in) :: &                  !
     243             :                        rat_h2oco2,rat_h2oco2_1, &
     244             :                        rat_h2oo3,rat_h2oo3_1, & !    Dimensions: (nlayers)
     245             :                        rat_h2on2o,rat_h2on2o_1, &
     246             :                        rat_h2och4,rat_h2och4_1, &
     247             :                        rat_n2oco2,rat_n2oco2_1, &
     248             :                        rat_o3co2,rat_o3co2_1
     249             : 
     250             : ! ----- Output -----
     251             :       real(kind=r8), intent(out) :: fracs(nlayers,ngptlw)        ! planck fractions
     252             :                                                       !    Dimensions: (nlayers,ngptlw)
     253             :       real(kind=r8), intent(out) :: taug(nlayers,ngptlw)         ! gaseous optical depth 
     254             :                                                       !    Dimensions: (nlayers,ngptlw)
     255             : 
     256             : ! Calculate gaseous optical depth and planck fractions for each spectral band.
     257             : 
     258      486000 :       call taugb1
     259      486000 :       call taugb2
     260      486000 :       call taugb3
     261      486000 :       call taugb4
     262      486000 :       call taugb5
     263      486000 :       call taugb6
     264      486000 :       call taugb7
     265      486000 :       call taugb8
     266      486000 :       call taugb9
     267      486000 :       call taugb10
     268      486000 :       call taugb11
     269      486000 :       call taugb12
     270      486000 :       call taugb13
     271      486000 :       call taugb14
     272      486000 :       call taugb15
     273      486000 :       call taugb16
     274             : 
     275             :       contains
     276             : 
     277             : !----------------------------------------------------------------------------
     278      486000 :       subroutine taugb1
     279             : !----------------------------------------------------------------------------
     280             : 
     281             : ! ------- Modifications -------
     282             : !  Written by Eli J. Mlawer, Atmospheric & Environmental Research.
     283             : !  Revised by Michael J. Iacono, Atmospheric & Environmental Research.
     284             : !
     285             : !     band 1:  10-350 cm-1 (low key - h2o; low minor - n2)
     286             : !                          (high key - h2o; high minor - n2)
     287             : !
     288             : !     note: previous versions of rrtm band 1: 
     289             : !           10-250 cm-1 (low - h2o; high - h2o)
     290             : !----------------------------------------------------------------------------
     291             : 
     292             : ! ------- Modules -------
     293             : 
     294             :       use parrrtm, only : ng1
     295             :       use rrlw_kg01, only : fracrefa, fracrefb, absa, absb, &
     296             :                             ka_mn2, kb_mn2, selfref, forref
     297             : 
     298             : ! ------- Declarations -------
     299             : 
     300             : ! Local 
     301      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     302      972000 :       real(kind=r8) ::  corradj(nlayers), scalen2(nlayers), tauself, taufor, taun2
     303             : 
     304             : 
     305             : ! Minor gas mapping levels:
     306             : !     lower - n2, p = 142.5490 mbar, t = 215.70 k
     307             : !     upper - n2, p = 142.5490 mbar, t = 215.70 k
     308             : 
     309             : ! Compute the optical depth by interpolating in ln(pressure) and 
     310             : ! temperature.  Below laytrop, the water vapor self-continuum and
     311             : ! foreign continuum is interpolated (in temperature) separately.
     312             : 
     313             : ! Lower atmosphere loop
     314    23763448 :       do lay = 1, laytrop
     315             : 
     316    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(1) + 1
     317    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(1) + 1
     318             :          !inds = indself(lay)
     319             :          !indf = indfor(lay)
     320             :          !indm = indminor(lay)
     321             :          !pp = pavel(lay)
     322    23277448 :          corradj(lay) =  1.
     323    23763448 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     324             :       enddo
     325    23763448 :       do lay = 1, laytrop
     326    23763448 :          if (pavel(lay) .lt. 250._r8) then
     327     6856962 :             corradj(lay) = 1._r8 - 0.15_r8 * (250._r8-pavel(lay)) / 154.4_r8
     328             :          endif
     329             :       enddo
     330    23763448 :       do lay = 1, laytrop
     331   256537928 :          do ig = 1, ng1
     332   698323440 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     333   931097920 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     334   465548960 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     335   465548960 :                  (forref(indfor(lay)+1,ig) -  forref(indfor(lay),ig))) 
     336   465548960 :             taun2 = scalen2(lay)*(ka_mn2(indminor(lay),ig) + & 
     337   465548960 :                  minorfrac(lay) * (ka_mn2(indminor(lay)+1,ig) - ka_mn2(indminor(lay),ig)))
     338   232774480 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     339   465548960 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     340   465548960 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     341   465548960 :                  fac01(lay) * absa(ind1(lay),ig) + &
     342   465548960 :                  fac11(lay) * absa(ind1(lay)+1,ig)) & 
     343   931097920 :                  + tauself + taufor + taun2)
     344   256051928 :              fracs(lay,ig) = fracrefa(ig)
     345             :          enddo
     346             :       enddo
     347             : 
     348             : ! Upper atmosphere loop
     349    22406552 :       do lay = laytrop+1, nlayers
     350             : 
     351    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(1) + 1
     352    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(1) + 1
     353             :          !indf = indfor(lay)
     354             :          !indm = indminor(lay)
     355             :          !pp = pavel(lay)
     356    21920552 :          corradj(lay) =  1._r8 - 0.15_r8 * (pavel(lay) / 95.6_r8)
     357             : 
     358    22406552 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     359             :       enddo
     360    22406552 :       do lay = laytrop+1, nlayers
     361   241612072 :          do ig = 1, ng1
     362   657616560 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     363   876822080 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     364   438411040 :             taun2 = scalen2(lay)*(kb_mn2(indminor(lay),ig) + & 
     365   438411040 :                  minorfrac(lay) * (kb_mn2(indminor(lay)+1,ig) - kb_mn2(indminor(lay),ig)))
     366   219205520 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     367   438411040 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     368   438411040 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     369   438411040 :                  fac01(lay) * absb(ind1(lay),ig) + &
     370   438411040 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &  
     371   876822080 :                  + taufor + taun2)
     372   241126072 :             fracs(lay,ig) = fracrefb(ig)
     373             :          enddo
     374             :       enddo
     375             : 
     376      486000 :       end subroutine taugb1
     377             : 
     378             : !----------------------------------------------------------------------------
     379      486000 :       subroutine taugb2
     380             : !----------------------------------------------------------------------------
     381             : !
     382             : !     band 2:  350-500 cm-1 (low key - h2o; high key - h2o)
     383             : !
     384             : !     note: previous version of rrtm band 2: 
     385             : !           250 - 500 cm-1 (low - h2o; high - h2o)
     386             : !----------------------------------------------------------------------------
     387             : 
     388             : ! ------- Modules -------
     389             : 
     390             :       use parrrtm, only : ng2, ngs1
     391             :       use rrlw_kg02, only : fracrefa, fracrefb, absa, absb, &
     392             :                             selfref, forref
     393             : 
     394             : ! ------- Declarations -------
     395             : 
     396             : ! Local 
     397      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     398      972000 :       real(kind=r8) :: corradj(nlayers), tauself, taufor
     399             : 
     400             : 
     401             : ! Compute the optical depth by interpolating in ln(pressure) and 
     402             : ! temperature.  Below laytrop, the water vapor self-continuum and
     403             : ! foreign continuum is interpolated (in temperature) separately.
     404             : 
     405             : ! Lower atmosphere loop
     406    23763448 :       do lay = 1, laytrop
     407             : 
     408    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(2) + 1
     409    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(2) + 1
     410             :          !inds = indself(lay)
     411             :          !indf = indfor(lay)
     412             :          !pp = pavel(lay)
     413    23763448 :          corradj(lay) = 1._r8 - .05_r8 * (pavel(lay) - 100._r8) / 900._r8
     414             :       enddo
     415    23763448 :       do lay = 1, laytrop
     416   303092824 :          do ig = 1, ng2
     417   837988128 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     418  1117317504 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     419   558658752 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     420   558658752 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     421   279329376 :             taug(lay,ngs1+ig) = corradj(lay) * (colh2o(lay) * &
     422   558658752 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     423   558658752 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     424   558658752 :                  fac01(lay) * absa(ind1(lay),ig) + &
     425   558658752 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
     426  1396646880 :                  + tauself + taufor)
     427   302606824 :             fracs(lay,ngs1+ig) = fracrefa(ig)
     428             :          enddo
     429             :       enddo
     430             : 
     431             : ! Upper atmosphere loop
     432    22406552 :       do lay = laytrop+1, nlayers
     433             : 
     434    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(2) + 1
     435    22406552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(2) + 1
     436             :          !indf = indfor(lay)
     437             :       enddo
     438    22406552 :       do lay = laytrop+1, nlayers
     439   285453176 :          do ig = 1, ng2
     440   789139872 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + &
     441  1052186496 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     442   263046624 :             taug(lay,ngs1+ig) = colh2o(lay) * &
     443   526093248 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     444   526093248 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     445   526093248 :                  fac01(lay) * absb(ind1(lay),ig) + &
     446   526093248 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
     447  1315233120 :                  + taufor
     448   284967176 :             fracs(lay,ngs1+ig) = fracrefb(ig)
     449             :          enddo
     450             :       enddo
     451             : 
     452      486000 :       end subroutine taugb2
     453             : 
     454             : !----------------------------------------------------------------------------
     455      486000 :       subroutine taugb3
     456             : !----------------------------------------------------------------------------
     457             : !
     458             : !     band 3:  500-630 cm-1 (low key - h2o,co2; low minor - n2o)
     459             : !                           (high key - h2o,co2; high minor - n2o)
     460             : !----------------------------------------------------------------------------
     461             : 
     462             : ! ------- Modules -------
     463             : 
     464             :       use parrrtm, only : ng3, ngs2
     465             :       use rrlw_ref, only : chi_mls
     466             :       use rrlw_kg03, only : fracrefa, fracrefb, absa, absb,&
     467             :                             ka_mn2o, kb_mn2o, selfref, forref
     468             : 
     469             : ! ------- Declarations -------
     470             : 
     471             : ! Local 
     472      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers),  ig
     473      972000 :       integer, dimension(nlayers) :: js, js1, jmn2o, jpl
     474      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
     475      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
     476      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2o, specparm_mn2o, specmult_mn2o, &
     477      972000 :                        fmn2o, fmn2omf, chi_n2o, ratn2o, adjfac, adjcoln2o
     478      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
     479             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
     480      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
     481      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
     482             :       real(kind=r8) :: tauself, taufor, n2om1, n2om2, absn2o
     483             :       real(kind=r8) :: refrat_planck_a, refrat_planck_b, refrat_m_a, refrat_m_b
     484             :       real(kind=r8) :: tau_major, tau_major1
     485             : 
     486             : 
     487             : ! Minor gas mapping levels:
     488             : !     lower - n2o, p = 706.272 mbar, t = 278.94 k
     489             : !     upper - n2o, p = 95.58 mbar, t = 215.7 k
     490             : 
     491             : !  P = 212.725 mb
     492      486000 :       refrat_planck_a = chi_mls(1,9)/chi_mls(2,9)
     493             : 
     494             : !  P = 95.58 mb
     495      486000 :       refrat_planck_b = chi_mls(1,13)/chi_mls(2,13)
     496             : 
     497             : !  P = 706.270mb
     498      486000 :       refrat_m_a = chi_mls(1,3)/chi_mls(2,3)
     499             : 
     500             : !  P = 95.58 mb 
     501      486000 :       refrat_m_b = chi_mls(1,13)/chi_mls(2,13)
     502             : 
     503             : ! Compute the optical depth by interpolating in ln(pressure) and 
     504             : ! temperature, and appropriate species.  Below laytrop, the water vapor 
     505             : ! self-continuum and foreign continuum is interpolated (in temperature) 
     506             : ! separately.
     507             : 
     508             : ! Lower atmosphere loop
     509    23763448 :       do lay = 1, laytrop
     510             : 
     511    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     512    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     513    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     514    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
     515    23277448 :          js(lay) = 1 + int(specmult(lay))
     516    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)        
     517             : 
     518    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     519    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     520    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     521    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
     522    23277448 :          js1(lay) = 1 + int(specmult1(lay))
     523    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     524             : 
     525    23277448 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
     526    23277448 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     527    23277448 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     528    23277448 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
     529    23277448 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     530    23277448 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     531    23277448 :          fmn2omf(lay) = minorfrac(lay)*fmn2o(lay)
     532             : !  In atmospheres where the amount of N2O is too great to be considered
     533             : !  a minor species, adjust the column amount of N2O by an empirical factor 
     534             : !  to obtain the proper contribution.
     535    23277448 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     536    23763448 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     537             :       enddo
     538    23763448 :       do lay = 1, laytrop
     539    23763448 :          if (ratn2o(lay) .gt. 1.5_r8) then
     540           0 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
     541           0 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
     542             :          else
     543    23277448 :             adjcoln2o(lay) = coln2o(lay)
     544             :          endif
     545             :       enddo
     546    23763448 :       do lay = 1, laytrop
     547             : 
     548    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     549    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     550    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     551    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     552    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
     553    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     554             : 
     555    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(3) + js(lay)
     556    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(3) + js1(lay)
     557             :       enddo
     558    23763448 :       do lay = 1, laytrop
     559             : 
     560    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
     561     3737204 :             p = fs(lay) - 1
     562     3737204 :             p4 = p**4
     563     3737204 :             fk0 = p4
     564     3737204 :             fk1 = 1 - p - 2.0_r8*p4
     565     3737204 :             fk2 = p + p4
     566     3737204 :             fac000(lay) = fk0*fac00(lay)
     567     3737204 :             fac100(lay) = fk1*fac00(lay)
     568     3737204 :             fac200(lay) = fk2*fac00(lay)
     569     3737204 :             fac010(lay) = fk0*fac10(lay)
     570     3737204 :             fac110(lay) = fk1*fac10(lay)
     571     3737204 :             fac210(lay) = fk2*fac10(lay)
     572    19540244 :          else if (specparm(lay) .gt. 0.875_r8) then
     573       31557 :             p = -fs(lay) 
     574       31557 :             p4 = p**4
     575       31557 :             fk0 = p4
     576       31557 :             fk1 = 1 - p - 2.0_r8*p4
     577       31557 :             fk2 = p + p4
     578       31557 :             fac000(lay) = fk0*fac00(lay)
     579       31557 :             fac100(lay) = fk1*fac00(lay)
     580       31557 :             fac200(lay) = fk2*fac00(lay)
     581       31557 :             fac010(lay) = fk0*fac10(lay)
     582       31557 :             fac110(lay) = fk1*fac10(lay)
     583       31557 :             fac210(lay) = fk2*fac10(lay)
     584             :          else
     585    19508687 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     586    19508687 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     587    19508687 :             fac100(lay) = fs(lay) * fac00(lay)
     588    19508687 :             fac110(lay) = fs(lay) * fac10(lay)
     589             :          endif
     590    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
     591     1487462 :             p = fs1(lay) - 1
     592     1487462 :             p4 = p**4
     593     1487462 :             fk0 = p4
     594     1487462 :             fk1 = 1 - p - 2.0_r8*p4
     595     1487462 :             fk2 = p + p4
     596     1487462 :             fac001(lay) = fk0*fac01(lay)
     597     1487462 :             fac101(lay) = fk1*fac01(lay)
     598     1487462 :             fac201(lay) = fk2*fac01(lay)
     599     1487462 :             fac011(lay) = fk0*fac11(lay)
     600     1487462 :             fac111(lay) = fk1*fac11(lay)
     601     1487462 :             fac211(lay) = fk2*fac11(lay)
     602    21789986 :          else if (specparm1(lay) .gt. 0.875_r8) then
     603      675028 :             p = -fs1(lay) 
     604      675028 :             p4 = p**4
     605      675028 :             fk0 = p4
     606      675028 :             fk1 = 1 - p - 2.0_r8*p4
     607      675028 :             fk2 = p + p4
     608      675028 :             fac001(lay) = fk0*fac01(lay)
     609      675028 :             fac101(lay) = fk1*fac01(lay)
     610      675028 :             fac201(lay) = fk2*fac01(lay)
     611      675028 :             fac011(lay) = fk0*fac11(lay)
     612      675028 :             fac111(lay) = fk1*fac11(lay)
     613      675028 :             fac211(lay) = fk2*fac11(lay)
     614             :          else
     615    21114958 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     616    21114958 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     617    21114958 :             fac101(lay) = fs1(lay) * fac01(lay)
     618    21114958 :             fac111(lay) = fs1(lay) * fac11(lay)
     619             :          endif
     620             :       enddo
     621    23763448 :       do lay = 1, laytrop
     622             : 
     623   396202616 :          do ig = 1, ng3
     624  1117317504 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     625  1489756672 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     626   744878336 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     627   744878336 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     628  1117317504 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     629  1117317504 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
     630   372439168 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     631   372439168 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     632   372439168 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     633             : 
     634   372439168 :             if (specparm(lay) .lt. 0.125_r8) then
     635             :                tau_major = speccomb(lay) * &
     636    59795264 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     637    59795264 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     638    59795264 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     639    59795264 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     640    59795264 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     641   358771584 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     642   312643904 :             else if (specparm(lay) .gt. 0.875_r8) then
     643             :                tau_major = speccomb(lay) * &
     644      504912 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     645      504912 :                     fac100(lay) * absa(ind0(lay),ig) + &
     646      504912 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     647      504912 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     648      504912 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     649     3029472 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     650             :             else
     651             :                tau_major = speccomb(lay) * &
     652   312138992 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     653   312138992 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     654   312138992 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     655  1248555968 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     656             :             endif
     657             : 
     658   372439168 :             if (specparm1(lay) .lt. 0.125_r8) then
     659             :                tau_major1 = speccomb1(lay) * &
     660    23799392 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     661    23799392 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     662    23799392 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     663    23799392 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     664    23799392 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     665   142796352 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     666   348639776 :             else if (specparm1(lay) .gt. 0.875_r8) then
     667             :                tau_major1 = speccomb1(lay) * &
     668    10800448 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     669    10800448 :                     fac101(lay) * absa(ind1(lay),ig) + &
     670    10800448 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     671    10800448 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     672    10800448 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     673    64802688 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     674             :             else
     675             :                tau_major1 = speccomb1(lay) * &
     676   337839328 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     677   337839328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     678   337839328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     679  1351357312 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     680             :             endif
     681             : 
     682   372439168 :             taug(lay,ngs2+ig) = tau_major + tau_major1 &
     683             :                  + tauself + taufor &
     684   372439168 :                  + adjcoln2o(lay)*absn2o
     685   744878336 :             fracs(lay,ngs2+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     686   768155784 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     687             :          enddo
     688             :       enddo
     689             : 
     690             : ! Upper atmosphere loop
     691    22406552 :       do lay = laytrop+1, nlayers
     692             : 
     693    21920552 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     694    21920552 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     695    21920552 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     696    21920552 :          specmult(lay) = 4._r8*(specparm(lay))
     697    21920552 :          js(lay) = 1 + int(specmult(lay))
     698    21920552 :          fs(lay) = mod(specmult(lay),1.0_r8)
     699             : 
     700    21920552 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     701    21920552 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     702    21920552 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     703    21920552 :          specmult1(lay) = 4._r8*(specparm1(lay))
     704    21920552 :          js1(lay) = 1 + int(specmult1(lay))
     705    21920552 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     706             : 
     707    21920552 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     708    21920552 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     709    21920552 :          fac100(lay) = fs(lay) * fac00(lay)
     710    21920552 :          fac110(lay) = fs(lay) * fac10(lay)
     711    21920552 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     712    21920552 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     713    21920552 :          fac101(lay) = fs1(lay) * fac01(lay)
     714    21920552 :          fac111(lay) = fs1(lay) * fac11(lay)
     715             : 
     716    21920552 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_b*colco2(lay)
     717    21920552 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     718    21920552 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     719    21920552 :          specmult_mn2o(lay) = 4._r8*specparm_mn2o(lay)
     720    21920552 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     721    21920552 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     722    21920552 :          fmn2omf(lay) = minorfrac(lay)*fmn2o(lay)
     723             : !  In atmospheres where the amount of N2O is too great to be considered
     724             : !  a minor species, adjust the column amount of N2O by an empirical factor 
     725             : !  to obtain the proper contribution.
     726    21920552 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     727    22406552 :          ratn2o(lay) = 1.e20*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     728             :       enddo
     729    22406552 :       do lay = laytrop+1, nlayers
     730    22406552 :          if (ratn2o(lay) .gt. 1.5_r8) then
     731    12971441 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
     732    12971441 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
     733             :          else
     734     8949111 :             adjcoln2o(lay) = coln2o(lay)
     735             :          endif
     736             :       enddo
     737    22406552 :       do lay = laytrop+1, nlayers
     738             : 
     739    21920552 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_b*colco2(lay)
     740    21920552 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     741    21920552 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     742    21920552 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
     743    21920552 :          jpl(lay)= 1 + int(specmult_planck(lay))
     744    21920552 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     745             : 
     746    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(3) + js(lay)
     747    22406552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(3) + js1(lay)
     748             :       enddo
     749    22406552 :       do lay = laytrop+1, nlayers
     750             : 
     751   373135384 :          do ig = 1, ng3
     752  1052186496 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     753  1402915328 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     754  1052186496 :             n2om1 = kb_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     755  1052186496 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay),ig)-kb_mn2o(jmn2o(lay),indminor(lay),ig))
     756   350728832 :             n2om2 = kb_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     757   350728832 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig)-kb_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     758   350728832 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     759   350728832 :             taug(lay,ngs2+ig) = speccomb(lay) * &
     760   350728832 :                 (fac000(lay) * absb(ind0(lay),ig) + &
     761   350728832 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
     762   350728832 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
     763   350728832 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
     764             :                 + speccomb1(lay) * &
     765   350728832 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
     766   350728832 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
     767   350728832 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
     768   350728832 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
     769             :                 + taufor &
     770  3156559488 :                 + adjcoln2o(lay)*absn2o
     771   701457664 :             fracs(lay,ngs2+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
     772   723378216 :                 (fracrefb(ig,jpl(lay)+1)-fracrefb(ig,jpl(lay)))
     773             :          enddo
     774             :       enddo
     775             : 
     776      486000 :       end subroutine taugb3
     777             : 
     778             : !----------------------------------------------------------------------------
     779      486000 :       subroutine taugb4
     780             : !----------------------------------------------------------------------------
     781             : !
     782             : !     band 4:  630-700 cm-1 (low key - h2o,co2; high key - o3,co2)
     783             : !----------------------------------------------------------------------------
     784             : 
     785             : ! ------- Modules -------
     786             : 
     787             :       use parrrtm, only : ng4, ngs3
     788             :       use rrlw_ref, only : chi_mls
     789             :       use rrlw_kg04, only : fracrefa, fracrefb, absa, absb, &
     790             :                             selfref, forref
     791             : 
     792             : ! ------- Declarations -------
     793             : 
     794             : ! Local 
     795      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     796      972000 :       integer, dimension(nlayers) :: js, js1, jpl
     797      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
     798      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
     799      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
     800             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
     801      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
     802      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
     803             :       real(kind=r8) :: tauself, taufor
     804             :       real(kind=r8) :: refrat_planck_a, refrat_planck_b
     805             :       real(kind=r8) :: tau_major, tau_major1
     806             : 
     807             : 
     808             : ! P =   142.5940 mb
     809      486000 :       refrat_planck_a = chi_mls(1,11)/chi_mls(2,11)
     810             : 
     811             : ! P = 95.58350 mb
     812      486000 :       refrat_planck_b = chi_mls(3,13)/chi_mls(2,13)
     813             : 
     814             : ! Compute the optical depth by interpolating in ln(pressure) and 
     815             : ! temperature, and appropriate species.  Below laytrop, the water 
     816             : ! vapor self-continuum and foreign continuum is interpolated (in temperature) 
     817             : ! separately.
     818             : 
     819             : ! Lower atmosphere loop
     820    23763448 :       do lay = 1, laytrop
     821             : 
     822    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     823    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     824    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     825    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
     826    23277448 :          js(lay) = 1 + int(specmult(lay))
     827    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
     828             : 
     829    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     830    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     831    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     832    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
     833    23277448 :          js1(lay) = 1 + int(specmult1(lay))
     834    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     835             : 
     836    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     837    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     838    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     839    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     840    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
     841    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     842             : 
     843    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(4) + js(lay)
     844    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(4) + js1(lay)
     845             :       enddo
     846    23763448 :       do lay = 1, laytrop
     847    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
     848     3737204 :             p = fs(lay) - 1
     849     3737204 :             p4 = p**4
     850     3737204 :             fk0 = p4
     851     3737204 :             fk1 = 1 - p - 2.0_r8*p4
     852     3737204 :             fk2 = p + p4
     853     3737204 :             fac000(lay) = fk0*fac00(lay)
     854     3737204 :             fac100(lay) = fk1*fac00(lay)
     855     3737204 :             fac200(lay) = fk2*fac00(lay)
     856     3737204 :             fac010(lay) = fk0*fac10(lay)
     857     3737204 :             fac110(lay) = fk1*fac10(lay)
     858     3737204 :             fac210(lay) = fk2*fac10(lay)
     859    19540244 :          else if (specparm(lay) .gt. 0.875_r8) then
     860       31557 :             p = -fs(lay) 
     861       31557 :             p4 = p**4
     862       31557 :             fk0 = p4
     863       31557 :             fk1 = 1 - p - 2.0_r8*p4
     864       31557 :             fk2 = p + p4
     865       31557 :             fac000(lay) = fk0*fac00(lay)
     866       31557 :             fac100(lay) = fk1*fac00(lay)
     867       31557 :             fac200(lay) = fk2*fac00(lay)
     868       31557 :             fac010(lay) = fk0*fac10(lay)
     869       31557 :             fac110(lay) = fk1*fac10(lay)
     870       31557 :             fac210(lay) = fk2*fac10(lay)
     871             :          else
     872    19508687 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     873    19508687 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     874    19508687 :             fac100(lay) = fs(lay) * fac00(lay)
     875    19508687 :             fac110(lay) = fs(lay) * fac10(lay)
     876             :          endif
     877             : 
     878    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
     879     1487462 :             p = fs1(lay) - 1
     880     1487462 :             p4 = p**4
     881     1487462 :             fk0 = p4
     882     1487462 :             fk1 = 1 - p - 2.0_r8*p4
     883     1487462 :             fk2 = p + p4
     884     1487462 :             fac001(lay) = fk0*fac01(lay)
     885     1487462 :             fac101(lay) = fk1*fac01(lay)
     886     1487462 :             fac201(lay) = fk2*fac01(lay)
     887     1487462 :             fac011(lay) = fk0*fac11(lay)
     888     1487462 :             fac111(lay) = fk1*fac11(lay)
     889     1487462 :             fac211(lay) = fk2*fac11(lay)
     890    21789986 :          else if (specparm1(lay) .gt. 0.875_r8) then
     891      675028 :             p = -fs1(lay) 
     892      675028 :             p4 = p**4
     893      675028 :             fk0 = p4
     894      675028 :             fk1 = 1 - p - 2.0_r8*p4
     895      675028 :             fk2 = p + p4
     896      675028 :             fac001(lay) = fk0*fac01(lay)
     897      675028 :             fac101(lay) = fk1*fac01(lay)
     898      675028 :             fac201(lay) = fk2*fac01(lay)
     899      675028 :             fac011(lay) = fk0*fac11(lay)
     900      675028 :             fac111(lay) = fk1*fac11(lay)
     901      675028 :             fac211(lay) = fk2*fac11(lay)
     902             :          else
     903    21114958 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     904    21114958 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     905    21114958 :             fac101(lay) = fs1(lay) * fac01(lay)
     906    21114958 :             fac111(lay) = fs1(lay) * fac11(lay)
     907             :          endif
     908             :       enddo
     909    23763448 :       do lay = 1, laytrop
     910             : 
     911   349647720 :          do ig = 1, ng4
     912   977652816 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     913  1303537088 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     914   651768544 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     915   651768544 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     916             : 
     917   325884272 :             if (specparm(lay) .lt. 0.125_r8) then
     918             :                tau_major = speccomb(lay) * &
     919    52320856 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     920    52320856 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     921    52320856 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     922    52320856 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     923    52320856 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     924   313925136 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     925   273563416 :             else if (specparm(lay) .gt. 0.875_r8) then
     926             :                tau_major = speccomb(lay) * &
     927      441798 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     928      441798 :                     fac100(lay) * absa(ind0(lay),ig) + &
     929      441798 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     930      441798 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     931      441798 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     932     2650788 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     933             :             else
     934             :                tau_major = speccomb(lay) * &
     935   273121618 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     936   273121618 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     937   273121618 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     938  1092486472 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     939             :             endif
     940             : 
     941   325884272 :             if (specparm1(lay) .lt. 0.125_r8) then
     942             :                tau_major1 = speccomb1(lay) * &
     943    20824468 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     944    20824468 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     945    20824468 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     946    20824468 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     947    20824468 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     948   124946808 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     949   305059804 :             else if (specparm1(lay) .gt. 0.875_r8) then
     950             :                tau_major1 = speccomb1(lay) * &
     951     9450392 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     952     9450392 :                     fac101(lay) * absa(ind1(lay),ig) + &
     953     9450392 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     954     9450392 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     955     9450392 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     956    56702352 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     957             :             else
     958             :                tau_major1 = speccomb1(lay) * &
     959   295609412 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     960   295609412 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     961   295609412 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     962  1182437648 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     963             :             endif
     964             : 
     965   325884272 :             taug(lay,ngs3+ig) = tau_major + tau_major1 &
     966   325884272 :                  + tauself + taufor
     967   651768544 :             fracs(lay,ngs3+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     968   675045992 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     969             :          enddo
     970             :       enddo
     971             : 
     972             : ! Upper atmosphere loop
     973    22406552 :       do lay = laytrop+1, nlayers
     974             : 
     975    21920552 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
     976    21920552 :          specparm(lay) = colo3(lay)/speccomb(lay)
     977    21920552 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     978    21920552 :          specmult(lay) = 4._r8*(specparm(lay))
     979    21920552 :          js(lay) = 1 + int(specmult(lay))
     980    21920552 :          fs(lay) = mod(specmult(lay),1.0_r8)
     981             : 
     982    21920552 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
     983    21920552 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
     984    21920552 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     985    21920552 :          specmult1(lay) = 4._r8*(specparm1(lay))
     986    21920552 :          js1(lay) = 1 + int(specmult1(lay))
     987    21920552 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     988             : 
     989    21920552 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     990    21920552 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     991    21920552 :          fac100(lay) = fs(lay) * fac00(lay)
     992    21920552 :          fac110(lay) = fs(lay) * fac10(lay)
     993    21920552 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     994    21920552 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     995    21920552 :          fac101(lay) = fs1(lay) * fac01(lay)
     996    21920552 :          fac111(lay) = fs1(lay) * fac11(lay)
     997             : 
     998    21920552 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
     999    21920552 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1000    21920552 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1001    21920552 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1002    21920552 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1003    21920552 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1004             : 
    1005    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(4) + js(lay)
    1006    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(4) + js1(lay)
    1007   328808280 :          do ig = 1, ng4
    1008   306887728 :             taug(lay,ngs3+ig) =  speccomb(lay) * &
    1009   613775456 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1010   306887728 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1011   306887728 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1012   306887728 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1013             :                 + speccomb1(lay) * &
    1014   306887728 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
    1015   306887728 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1016   306887728 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1017  3068877280 :                 fac111(lay) * absb(ind1(lay)+6,ig))
    1018   613775456 :             fracs(lay,ngs3+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1019   635696008 :                 (fracrefb(ig,jpl(lay)+1)-fracrefb(ig,jpl(lay)))
    1020             :          enddo
    1021             : 
    1022             : ! Empirical modification to code to improve stratospheric cooling rates
    1023             : ! for co2.  Revised to apply weighting for g-point reduction in this band.
    1024             : 
    1025    21920552 :          taug(lay,ngs3+8)=taug(lay,ngs3+8)*0.92
    1026    21920552 :          taug(lay,ngs3+9)=taug(lay,ngs3+9)*0.88
    1027    21920552 :          taug(lay,ngs3+10)=taug(lay,ngs3+10)*1.07
    1028    21920552 :          taug(lay,ngs3+11)=taug(lay,ngs3+11)*1.1
    1029    21920552 :          taug(lay,ngs3+12)=taug(lay,ngs3+12)*0.99
    1030    21920552 :          taug(lay,ngs3+13)=taug(lay,ngs3+13)*0.88
    1031    22406552 :          taug(lay,ngs3+14)=taug(lay,ngs3+14)*0.943
    1032             : 
    1033             :       enddo
    1034             : 
    1035      486000 :       end subroutine taugb4
    1036             : 
    1037             : !----------------------------------------------------------------------------
    1038      486000 :       subroutine taugb5
    1039             : !----------------------------------------------------------------------------
    1040             : !
    1041             : !     band 5:  700-820 cm-1 (low key - h2o,co2; low minor - o3, ccl4)
    1042             : !                           (high key - o3,co2)
    1043             : !----------------------------------------------------------------------------
    1044             : 
    1045             : ! ------- Modules -------
    1046             : 
    1047             :       use parrrtm, only : ng5, ngs4
    1048             :       use rrlw_ref, only : chi_mls
    1049             :       use rrlw_kg05, only : fracrefa, fracrefb, absa, absb, &
    1050             :                             ka_mo3, selfref, forref, ccl4
    1051             : 
    1052             : ! ------- Declarations -------
    1053             : 
    1054             : ! Local 
    1055      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers),  ig
    1056      972000 :       integer, dimension(nlayers) :: js, js1, jmo3, jpl
    1057      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1058      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1059      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mo3, specparm_mo3, specmult_mo3, fmo3
    1060      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1061             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1062      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1063      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    1064             :       real(kind=r8) :: tauself, taufor, o3m1, o3m2, abso3
    1065             :       real(kind=r8) :: refrat_planck_a, refrat_planck_b, refrat_m_a
    1066             :       real(kind=r8) :: tau_major, tau_major1
    1067             : 
    1068             : 
    1069             : ! Minor gas mapping level :
    1070             : !     lower - o3, p = 317.34 mbar, t = 240.77 k
    1071             : !     lower - ccl4
    1072             : 
    1073             : ! Calculate reference ratio to be used in calculation of Planck
    1074             : ! fraction in lower/upper atmosphere.
    1075             : 
    1076             : ! P = 473.420 mb
    1077      486000 :       refrat_planck_a = chi_mls(1,5)/chi_mls(2,5)
    1078             : 
    1079             : ! P = 0.2369 mb
    1080      486000 :       refrat_planck_b = chi_mls(3,43)/chi_mls(2,43)
    1081             : 
    1082             : ! P = 317.3480
    1083      486000 :       refrat_m_a = chi_mls(1,7)/chi_mls(2,7)
    1084             : 
    1085             : ! Compute the optical depth by interpolating in ln(pressure) and 
    1086             : ! temperature, and appropriate species.  Below laytrop, the 
    1087             : ! water vapor self-continuum and foreign continuum is 
    1088             : ! interpolated (in temperature) separately.
    1089             : 
    1090             : ! Lower atmosphere loop
    1091    23763448 :       do lay = 1, laytrop
    1092             : 
    1093    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    1094    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1095    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1096    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    1097    23277448 :          js(lay) = 1 + int(specmult(lay))
    1098    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1099             : 
    1100    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    1101    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1102    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1103    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1104    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    1105    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1106             : 
    1107    23277448 :          speccomb_mo3(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
    1108    23277448 :          specparm_mo3(lay) = colh2o(lay)/speccomb_mo3(lay)
    1109    23277448 :          if (specparm_mo3(lay) .ge. oneminus) specparm_mo3(lay) = oneminus
    1110    23277448 :          specmult_mo3(lay) = 8._r8*specparm_mo3(lay)
    1111    23277448 :          jmo3(lay) = 1 + int(specmult_mo3(lay))
    1112    23277448 :          fmo3(lay) = mod(specmult_mo3(lay),1.0_r8)
    1113             : 
    1114    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    1115    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1116    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1117    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1118    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1119    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1120             : 
    1121    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(5) + js(lay)
    1122    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(5) + js1(lay)
    1123             :       enddo
    1124    23763448 :       do lay = 1, laytrop
    1125    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
    1126     3737204 :             p = fs(lay) - 1
    1127     3737204 :             p4 = p**4
    1128     3737204 :             fk0 = p4
    1129     3737204 :             fk1 = 1 - p - 2.0_r8*p4
    1130     3737204 :             fk2 = p + p4
    1131     3737204 :             fac000(lay) = fk0*fac00(lay)
    1132     3737204 :             fac100(lay) = fk1*fac00(lay)
    1133     3737204 :             fac200(lay) = fk2*fac00(lay)
    1134     3737204 :             fac010(lay) = fk0*fac10(lay)
    1135     3737204 :             fac110(lay) = fk1*fac10(lay)
    1136     3737204 :             fac210(lay) = fk2*fac10(lay)
    1137    19540244 :          else if (specparm(lay) .gt. 0.875_r8) then
    1138       31557 :             p = -fs(lay) 
    1139       31557 :             p4 = p**4
    1140       31557 :             fk0 = p4
    1141       31557 :             fk1 = 1 - p - 2.0_r8*p4
    1142       31557 :             fk2 = p + p4
    1143       31557 :             fac000(lay) = fk0*fac00(lay)
    1144       31557 :             fac100(lay) = fk1*fac00(lay)
    1145       31557 :             fac200(lay) = fk2*fac00(lay)
    1146       31557 :             fac010(lay) = fk0*fac10(lay)
    1147       31557 :             fac110(lay) = fk1*fac10(lay)
    1148       31557 :             fac210(lay) = fk2*fac10(lay)
    1149             :          else
    1150    19508687 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1151    19508687 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1152    19508687 :             fac100(lay) = fs(lay) * fac00(lay)
    1153    19508687 :             fac110(lay) = fs(lay) * fac10(lay)
    1154             :          endif
    1155             : 
    1156    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
    1157     1487462 :             p = fs1(lay) - 1
    1158     1487462 :             p4 = p**4
    1159     1487462 :             fk0 = p4
    1160     1487462 :             fk1 = 1 - p - 2.0_r8*p4
    1161     1487462 :             fk2 = p + p4
    1162     1487462 :             fac001(lay) = fk0*fac01(lay)
    1163     1487462 :             fac101(lay) = fk1*fac01(lay)
    1164     1487462 :             fac201(lay) = fk2*fac01(lay)
    1165     1487462 :             fac011(lay) = fk0*fac11(lay)
    1166     1487462 :             fac111(lay) = fk1*fac11(lay)
    1167     1487462 :             fac211(lay) = fk2*fac11(lay)
    1168    21789986 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1169      675028 :             p = -fs1(lay) 
    1170      675028 :             p4 = p**4
    1171      675028 :             fk0 = p4
    1172      675028 :             fk1 = 1 - p - 2.0_r8*p4
    1173      675028 :             fk2 = p + p4
    1174      675028 :             fac001(lay) = fk0*fac01(lay)
    1175      675028 :             fac101(lay) = fk1*fac01(lay)
    1176      675028 :             fac201(lay) = fk2*fac01(lay)
    1177      675028 :             fac011(lay) = fk0*fac11(lay)
    1178      675028 :             fac111(lay) = fk1*fac11(lay)
    1179      675028 :             fac211(lay) = fk2*fac11(lay)
    1180             :          else
    1181    21114958 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1182    21114958 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1183    21114958 :             fac101(lay) = fs1(lay) * fac01(lay)
    1184    21114958 :             fac111(lay) = fs1(lay) * fac11(lay)
    1185             :          endif
    1186             :       enddo
    1187    23763448 :       do lay = 1, laytrop
    1188             : 
    1189   396202616 :          do ig = 1, ng5
    1190  1117317504 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1191  1489756672 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1192   744878336 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1193   744878336 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1194  1117317504 :             o3m1 = ka_mo3(jmo3(lay),indminor(lay),ig) + fmo3(lay) * &
    1195  1117317504 :                  (ka_mo3(jmo3(lay)+1,indminor(lay),ig)-ka_mo3(jmo3(lay),indminor(lay),ig))
    1196   372439168 :             o3m2 = ka_mo3(jmo3(lay),indminor(lay)+1,ig) + fmo3(lay) * &
    1197   372439168 :                  (ka_mo3(jmo3(lay)+1,indminor(lay)+1,ig)-ka_mo3(jmo3(lay),indminor(lay)+1,ig))
    1198   372439168 :             abso3 = o3m1 + minorfrac(lay)*(o3m2-o3m1)
    1199             : 
    1200   372439168 :             if (specparm(lay) .lt. 0.125_r8) then
    1201             :                tau_major = speccomb(lay) * &
    1202    59795264 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1203    59795264 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1204    59795264 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1205    59795264 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1206    59795264 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1207   358771584 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1208   312643904 :             else if (specparm(lay) .gt. 0.875_r8) then
    1209             :                tau_major = speccomb(lay) * &
    1210      504912 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1211      504912 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1212      504912 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1213      504912 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1214      504912 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1215     3029472 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1216             :             else
    1217             :                tau_major = speccomb(lay) * &
    1218   312138992 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1219   312138992 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1220   312138992 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1221  1248555968 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1222             :             endif
    1223             : 
    1224   372439168 :             if (specparm1(lay) .lt. 0.125_r8) then
    1225             :                tau_major1 = speccomb1(lay) * &
    1226    23799392 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1227    23799392 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1228    23799392 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1229    23799392 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1230    23799392 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1231   142796352 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1232   348639776 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1233             :                tau_major1 = speccomb1(lay) * & 
    1234    10800448 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1235    10800448 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1236    10800448 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1237    10800448 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1238    10800448 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1239    64802688 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1240             :             else
    1241             :                tau_major1 = speccomb1(lay) * &
    1242   337839328 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1243   337839328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1244   337839328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1245  1351357312 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1246             :             endif
    1247             : 
    1248   372439168 :             taug(lay,ngs4+ig) = tau_major + tau_major1 &
    1249             :                  + tauself + taufor &
    1250   372439168 :                  + abso3*colo3(lay) &
    1251   744878336 :                  + wx(1,lay) * ccl4(ig)
    1252   744878336 :             fracs(lay,ngs4+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1253   768155784 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1254             :          enddo
    1255             :       enddo
    1256             : 
    1257             : ! Upper atmosphere loop
    1258    22406552 :       do lay = laytrop+1, nlayers
    1259             : 
    1260    21920552 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
    1261    21920552 :          specparm(lay) = colo3(lay)/speccomb(lay)
    1262    21920552 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1263    21920552 :          specmult(lay) = 4._r8*(specparm(lay))
    1264    21920552 :          js(lay) = 1 + int(specmult(lay))
    1265    21920552 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1266             : 
    1267    21920552 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
    1268    21920552 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
    1269    21920552 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1270    21920552 :          specmult1(lay) = 4._r8*(specparm1(lay))
    1271    21920552 :          js1(lay) = 1 + int(specmult1(lay))
    1272    21920552 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1273             : 
    1274    21920552 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1275    21920552 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1276    21920552 :          fac100(lay) = fs(lay) * fac00(lay)
    1277    21920552 :          fac110(lay) = fs(lay) * fac10(lay)
    1278    21920552 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1279    21920552 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1280    21920552 :          fac101(lay) = fs1(lay) * fac01(lay)
    1281    21920552 :          fac111(lay) = fs1(lay) * fac11(lay)
    1282             : 
    1283    21920552 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
    1284    21920552 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1285    21920552 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1286    21920552 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1287    21920552 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1288    21920552 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1289             : 
    1290    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(5) + js(lay)
    1291    22406552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(5) + js1(lay)
    1292             :       enddo
    1293    22406552 :       do lay = laytrop+1, nlayers
    1294   373135384 :          do ig = 1, ng5
    1295   701457664 :             taug(lay,ngs4+ig) = speccomb(lay) * &
    1296   701457664 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1297   350728832 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1298   350728832 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1299   350728832 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1300             :                 + speccomb1(lay) * &
    1301   350728832 :                 (fac001(lay) * absb(ind1(lay),ig) + &
    1302   350728832 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1303   350728832 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1304   350728832 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
    1305  4208745984 :                 + wx(1,lay) * ccl4(ig)
    1306   701457664 :             fracs(lay,ngs4+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1307   723378216 :                 (fracrefb(ig,jpl(lay)+1)-fracrefb(ig,jpl(lay)))
    1308             :          enddo
    1309             :       enddo
    1310             : 
    1311      486000 :       end subroutine taugb5
    1312             : 
    1313             : !----------------------------------------------------------------------------
    1314      486000 :       subroutine taugb6
    1315             : !----------------------------------------------------------------------------
    1316             : !
    1317             : !     band 6:  820-980 cm-1 (low key - h2o; low minor - co2)
    1318             : !                           (high key - nothing; high minor - cfc11, cfc12)
    1319             : !----------------------------------------------------------------------------
    1320             : 
    1321             : ! ------- Modules -------
    1322             : 
    1323             :       use parrrtm, only : ng6, ngs5
    1324             :       use rrlw_ref, only : chi_mls
    1325             :       use rrlw_kg06, only : fracrefa, absa, ka_mco2, &
    1326             :                             selfref, forref, cfc11adj, cfc12
    1327             : 
    1328             : ! ------- Declarations -------
    1329             : 
    1330             : ! Local 
    1331      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1332      972000 :       real(kind=r8), dimension(nlayers) :: chi_co2, ratco2, adjfac, adjcolco2
    1333             :       real(kind=r8) :: tauself, taufor, absco2
    1334             : 
    1335             : 
    1336             : ! Minor gas mapping level:
    1337             : !     lower - co2, p = 706.2720 mb, t = 294.2 k
    1338             : !     upper - cfc11, cfc12
    1339             : 
    1340             : ! Compute the optical depth by interpolating in ln(pressure) and
    1341             : ! temperature. The water vapor self-continuum and foreign continuum
    1342             : ! is interpolated (in temperature) separately.  
    1343             : 
    1344             : ! Lower atmosphere loop
    1345    23763448 :       do lay = 1, laytrop
    1346             : 
    1347             : ! In atmospheres where the amount of CO2 is too great to be considered
    1348             : ! a minor species, adjust the column amount of CO2 by an empirical factor 
    1349             : ! to obtain the proper contribution.
    1350    23277448 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1351    23763448 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1352             :       enddo
    1353    23763448 :       do lay = 1, laytrop
    1354    23763448 :          if (ratco2(lay) .gt. 3.0_r8) then
    1355           0 :             adjfac(lay) = 2.0_r8+(ratco2(lay)-2.0_r8)**0.77_r8
    1356           0 :             adjcolco2(lay) = adjfac(lay)*chi_mls(2,jp(lay)+1)*coldry(lay)*1.e-20_r8
    1357             :          else
    1358    23277448 :             adjcolco2(lay) = colco2(lay)
    1359             :          endif
    1360             :       enddo
    1361    23763448 :       do lay = 1, laytrop
    1362             : 
    1363    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(6) + 1
    1364    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(6) + 1
    1365             :       enddo
    1366    23763448 :       do lay = 1, laytrop
    1367             : 
    1368   209983032 :          do ig = 1, ng6
    1369   558658752 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1370   744878336 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1371   372439168 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1372   372439168 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1373   372439168 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1374   372439168 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1375   186219584 :             taug(lay,ngs5+ig) = colh2o(lay) * &
    1376   372439168 :                 (fac00(lay) * absa(ind0(lay),ig) + &
    1377   372439168 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1378   372439168 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1379   372439168 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    1380             :                  + tauself + taufor &
    1381             :                  + adjcolco2(lay) * absco2 &
    1382   186219584 :                  + wx(2,lay) * cfc11adj(ig) &
    1383   931097920 :                  + wx(3,lay) * cfc12(ig)
    1384   209497032 :             fracs(lay,ngs5+ig) = fracrefa(ig)
    1385             :          enddo
    1386             :       enddo
    1387             : 
    1388             : ! Upper atmosphere loop
    1389             : ! Nothing important goes on above laytrop in this band.
    1390    22406552 :       do lay = laytrop+1, nlayers
    1391             : 
    1392   197770968 :          do ig = 1, ng6
    1393   350728832 :             taug(lay,ngs5+ig) = 0.0_r8 &
    1394   350728832 :                  + wx(2,lay) * cfc11adj(ig) &
    1395   526093248 :                  + wx(3,lay) * cfc12(ig)
    1396   197284968 :             fracs(lay,ngs5+ig) = fracrefa(ig)
    1397             :          enddo
    1398             :       enddo
    1399             : 
    1400      486000 :       end subroutine taugb6
    1401             : 
    1402             : !----------------------------------------------------------------------------
    1403      486000 :       subroutine taugb7
    1404             : !----------------------------------------------------------------------------
    1405             : !
    1406             : !     band 7:  980-1080 cm-1 (low key - h2o,o3; low minor - co2)
    1407             : !                            (high key - o3; high minor - co2)
    1408             : !----------------------------------------------------------------------------
    1409             : 
    1410             : ! ------- Modules -------
    1411             : 
    1412             :       use parrrtm, only : ng7, ngs6
    1413             :       use rrlw_ref, only : chi_mls
    1414             :       use rrlw_kg07, only : fracrefa, fracrefb, absa, absb, &
    1415             :                             ka_mco2, kb_mco2, selfref, forref
    1416             : 
    1417             : ! ------- Declarations -------
    1418             : 
    1419             : ! Local 
    1420      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1421      972000 :       integer, dimension(nlayers) :: js, js1, jmco2, jpl
    1422      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1423      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1424      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mco2, specparm_mco2, specmult_mco2, fmco2
    1425      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1426             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1427      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1428      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    1429             :       real(kind=r8) :: tauself, taufor, co2m1, co2m2, absco2
    1430      972000 :       real(kind=r8), dimension(nlayers) :: chi_co2, ratco2, adjfac, adjcolco2
    1431             :       real(kind=r8) :: refrat_planck_a, refrat_m_a
    1432             :       real(kind=r8) :: tau_major, tau_major1
    1433             : 
    1434             : 
    1435             : ! Minor gas mapping level :
    1436             : !     lower - co2, p = 706.2620 mbar, t= 278.94 k
    1437             : !     upper - co2, p = 12.9350 mbar, t = 234.01 k
    1438             : 
    1439             : ! Calculate reference ratio to be used in calculation of Planck
    1440             : ! fraction in lower atmosphere.
    1441             : 
    1442             : ! P = 706.2620 mb
    1443      486000 :       refrat_planck_a = chi_mls(1,3)/chi_mls(3,3)
    1444             : 
    1445             : ! P = 706.2720 mb
    1446      486000 :       refrat_m_a = chi_mls(1,3)/chi_mls(3,3)
    1447             : 
    1448             : ! Compute the optical depth by interpolating in ln(pressure), 
    1449             : ! temperature, and appropriate species.  Below laytrop, the water
    1450             : ! vapor self-continuum and foreign continuum is interpolated 
    1451             : ! (in temperature) separately. 
    1452             : 
    1453             : ! Lower atmosphere loop
    1454    23763448 :       do lay = 1, laytrop
    1455             : 
    1456    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2oo3(lay)*colo3(lay)
    1457    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1458    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1459    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    1460    23277448 :          js(lay) = 1 + int(specmult(lay))
    1461    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1462             : 
    1463    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2oo3_1(lay)*colo3(lay)
    1464    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1465    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1466    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1467    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    1468    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1469             : 
    1470    23277448 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*colo3(lay)
    1471    23277448 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    1472    23277448 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    1473    23277448 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    1474             : 
    1475    23277448 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    1476    23277448 :          fmco2(lay) = mod(specmult_mco2(lay),1.0_r8)
    1477             : 
    1478             : !  In atmospheres where the amount of CO2 is too great to be considered
    1479             : !  a minor species, adjust the column amount of CO2 by an empirical factor 
    1480             : !  to obtain the proper contribution.
    1481    23277448 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1482    23277448 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1483    23277448 :          if (ratco2(lay) .gt. 3.0_r8) then
    1484           0 :             adjfac(lay) = 3.0_r8+(ratco2(lay)-3.0_r8)**0.79_r8
    1485           0 :             adjcolco2(lay) = adjfac(lay)*chi_mls(2,jp(lay)+1)*coldry(lay)*1.e-20_r8
    1486             :          else
    1487    23277448 :             adjcolco2(lay) = colco2(lay)
    1488             :          endif
    1489             : 
    1490    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colo3(lay)
    1491    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1492    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1493    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1494    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1495    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1496             : 
    1497    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(7) + js(lay)
    1498    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(7) + js1(lay)
    1499             : 
    1500    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
    1501     3446017 :             p = fs(lay) - 1
    1502     3446017 :             p4 = p**4
    1503     3446017 :             fk0 = p4
    1504     3446017 :             fk1 = 1 - p - 2.0_r8*p4
    1505     3446017 :             fk2 = p + p4
    1506     3446017 :             fac000(lay) = fk0*fac00(lay)
    1507     3446017 :             fac100(lay) = fk1*fac00(lay)
    1508     3446017 :             fac200(lay) = fk2*fac00(lay)
    1509     3446017 :             fac010(lay) = fk0*fac10(lay)
    1510     3446017 :             fac110(lay) = fk1*fac10(lay)
    1511     3446017 :             fac210(lay) = fk2*fac10(lay)
    1512    19831431 :          else if (specparm(lay) .gt. 0.875_r8) then
    1513     1893189 :             p = -fs(lay) 
    1514     1893189 :             p4 = p**4
    1515     1893189 :             fk0 = p4
    1516     1893189 :             fk1 = 1 - p - 2.0_r8*p4
    1517     1893189 :             fk2 = p + p4
    1518     1893189 :             fac000(lay) = fk0*fac00(lay)
    1519     1893189 :             fac100(lay) = fk1*fac00(lay)
    1520     1893189 :             fac200(lay) = fk2*fac00(lay)
    1521     1893189 :             fac010(lay) = fk0*fac10(lay)
    1522     1893189 :             fac110(lay) = fk1*fac10(lay)
    1523     1893189 :             fac210(lay) = fk2*fac10(lay)
    1524             :          else
    1525    17938242 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1526    17938242 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1527    17938242 :             fac100(lay) = fs(lay) * fac00(lay)
    1528    17938242 :             fac110(lay) = fs(lay) * fac10(lay)
    1529             :          endif
    1530    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
    1531     3446017 :             p = fs1(lay) - 1
    1532     3446017 :             p4 = p**4
    1533     3446017 :             fk0 = p4
    1534     3446017 :             fk1 = 1 - p - 2.0_r8*p4
    1535     3446017 :             fk2 = p + p4
    1536     3446017 :             fac001(lay) = fk0*fac01(lay)
    1537     3446017 :             fac101(lay) = fk1*fac01(lay)
    1538     3446017 :             fac201(lay) = fk2*fac01(lay)
    1539     3446017 :             fac011(lay) = fk0*fac11(lay)
    1540     3446017 :             fac111(lay) = fk1*fac11(lay)
    1541     3446017 :             fac211(lay) = fk2*fac11(lay)
    1542    19831431 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1543     4063738 :             p = -fs1(lay) 
    1544     4063738 :             p4 = p**4
    1545     4063738 :             fk0 = p4
    1546     4063738 :             fk1 = 1 - p - 2.0_r8*p4
    1547     4063738 :             fk2 = p + p4
    1548     4063738 :             fac001(lay) = fk0*fac01(lay)
    1549     4063738 :             fac101(lay) = fk1*fac01(lay)
    1550     4063738 :             fac201(lay) = fk2*fac01(lay)
    1551     4063738 :             fac011(lay) = fk0*fac11(lay)
    1552     4063738 :             fac111(lay) = fk1*fac11(lay)
    1553     4063738 :             fac211(lay) = fk2*fac11(lay)
    1554             :          else
    1555    15767693 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1556    15767693 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1557    15767693 :             fac101(lay) = fs1(lay) * fac01(lay)
    1558    15767693 :             fac111(lay) = fs1(lay) * fac11(lay)
    1559             :          endif
    1560             : 
    1561   303092824 :          do ig = 1, ng7
    1562   837988128 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1563   837988128 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1564   558658752 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1565   558658752 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1566   837988128 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    1567   837988128 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    1568   279329376 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    1569   279329376 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    1570   279329376 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    1571             : 
    1572   279329376 :             if (specparm(lay) .lt. 0.125_r8) then
    1573             :                tau_major = speccomb(lay) * &
    1574    41352204 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1575    41352204 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1576    41352204 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1577    41352204 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1578    41352204 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1579   248113224 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1580   237977172 :             else if (specparm(lay) .gt. 0.875_r8) then
    1581             :                tau_major = speccomb(lay) * &
    1582    22718268 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1583    22718268 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1584    22718268 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1585    22718268 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1586    22718268 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1587   136309608 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1588             :             else
    1589             :                tau_major = speccomb(lay) * &
    1590   215258904 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1591   215258904 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1592   215258904 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1593   861035616 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1594             :             endif
    1595             : 
    1596   279329376 :             if (specparm1(lay) .lt. 0.125_r8) then
    1597             :                tau_major1 = speccomb1(lay) * &
    1598    13103076 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1599    13103076 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1600    13103076 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1601    13103076 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1602    13103076 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1603    78618456 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1604   266226300 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1605             :                tau_major1 = speccomb1(lay) * &
    1606    48764856 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1607    48764856 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1608    48764856 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1609    48764856 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1610    48764856 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1611   292589136 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1612             :             else
    1613             :                tau_major1 = speccomb1(lay) * &
    1614   217461444 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
    1615   217461444 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1616   217461444 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1617   869845776 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1618             :             endif
    1619             : 
    1620   279329376 :             taug(lay,ngs6+ig) = tau_major + tau_major1 &
    1621             :                  + tauself + taufor &
    1622   279329376 :                  + adjcolco2(lay)*absco2
    1623   558658752 :             fracs(lay,ngs6+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1624   581936200 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1625             :          enddo
    1626             :       enddo
    1627             : 
    1628             : ! Upper atmosphere loop
    1629    22406552 :       do lay = laytrop+1, nlayers
    1630             : 
    1631             : !  In atmospheres where the amount of CO2 is too great to be considered
    1632             : !  a minor species, adjust the column amount of CO2 by an empirical factor 
    1633             : !  to obtain the proper contribution.
    1634    21920552 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1635    21920552 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1636    21920552 :          if (ratco2(lay) .gt. 3.0_r8) then
    1637           0 :             adjfac(lay) = 2.0_r8+(ratco2(lay)-2.0_r8)**0.79_r8
    1638           0 :             adjcolco2(lay) = adjfac(lay)*chi_mls(2,jp(lay)+1)*coldry(lay)*1.e-20_r8
    1639             :          else
    1640    21920552 :             adjcolco2(lay) = colco2(lay)
    1641             :          endif
    1642             : 
    1643    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(7) + 1
    1644    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(7) + 1
    1645             : 
    1646   284967176 :          do ig = 1, ng7
    1647   526093248 :             absco2 = kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1648   789139872 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig))
    1649   263046624 :             taug(lay,ngs6+ig) = colo3(lay) * &
    1650   526093248 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1651   526093248 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1652   526093248 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1653   526093248 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1654  1315233120 :                  + adjcolco2(lay) * absco2
    1655   284967176 :             fracs(lay,ngs6+ig) = fracrefb(ig)
    1656             :          enddo
    1657             : 
    1658             : ! Empirical modification to code to improve stratospheric cooling rates
    1659             : ! for o3.  Revised to apply weighting for g-point reduction in this band.
    1660             : 
    1661    21920552 :          taug(lay,ngs6+6)=taug(lay,ngs6+6)*0.92_r8
    1662    21920552 :          taug(lay,ngs6+7)=taug(lay,ngs6+7)*0.88_r8
    1663    21920552 :          taug(lay,ngs6+8)=taug(lay,ngs6+8)*1.07_r8
    1664    21920552 :          taug(lay,ngs6+9)=taug(lay,ngs6+9)*1.1_r8
    1665    21920552 :          taug(lay,ngs6+10)=taug(lay,ngs6+10)*0.99_r8
    1666    22406552 :          taug(lay,ngs6+11)=taug(lay,ngs6+11)*0.855_r8
    1667             : 
    1668             :       enddo
    1669             : 
    1670      486000 :       end subroutine taugb7
    1671             : 
    1672             : !----------------------------------------------------------------------------
    1673      486000 :       subroutine taugb8
    1674             : !----------------------------------------------------------------------------
    1675             : !
    1676             : !     band 8:  1080-1180 cm-1 (low key - h2o; low minor - co2,o3,n2o)
    1677             : !                             (high key - o3; high minor - co2, n2o)
    1678             : !----------------------------------------------------------------------------
    1679             : 
    1680             : ! ------- Modules -------
    1681             : 
    1682             :       use parrrtm, only : ng8, ngs7
    1683             :       use rrlw_ref, only : chi_mls
    1684             :       use rrlw_kg08, only : fracrefa, fracrefb, absa, absb, &
    1685             :                             ka_mco2, ka_mn2o, ka_mo3, kb_mco2, kb_mn2o, &
    1686             :                             selfref, forref, cfc12, cfc22adj
    1687             : 
    1688             : ! ------- Declarations -------
    1689             : 
    1690             : ! Local 
    1691      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1692             :       real(kind=r8) :: tauself, taufor, absco2, abso3, absn2o
    1693      972000 :       real(kind=r8), dimension(nlayers) :: chi_co2, ratco2, adjfac, adjcolco2
    1694             : 
    1695             : 
    1696             : ! Minor gas mapping level:
    1697             : !     lower - co2, p = 1053.63 mb, t = 294.2 k
    1698             : !     lower - o3,  p = 317.348 mb, t = 240.77 k
    1699             : !     lower - n2o, p = 706.2720 mb, t= 278.94 k
    1700             : !     lower - cfc12,cfc11
    1701             : !     upper - co2, p = 35.1632 mb, t = 223.28 k
    1702             : !     upper - n2o, p = 8.716e-2 mb, t = 226.03 k
    1703             : 
    1704             : ! Compute the optical depth by interpolating in ln(pressure) and 
    1705             : ! temperature, and appropriate species.  Below laytrop, the water vapor 
    1706             : ! self-continuum and foreign continuum is interpolated (in temperature) 
    1707             : ! separately.
    1708             : 
    1709             : ! Lower atmosphere loop
    1710    23763448 :       do lay = 1, laytrop
    1711             : 
    1712             : !  In atmospheres where the amount of CO2 is too great to be considered
    1713             : !  a minor species, adjust the column amount of CO2 by an empirical factor 
    1714             : !  to obtain the proper contribution.
    1715    23277448 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1716    23277448 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1717    23277448 :          if (ratco2(lay) .gt. 3.0_r8) then
    1718           0 :             adjfac(lay) = 2.0_r8+(ratco2(lay)-2.0_r8)**0.65_r8
    1719           0 :             adjcolco2(lay) = adjfac(lay)*chi_mls(2,jp(lay)+1)*coldry(lay)*1.e-20_r8
    1720             :          else
    1721    23277448 :             adjcolco2(lay) = colco2(lay)
    1722             :          endif
    1723             : 
    1724    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(8) + 1
    1725    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(8) + 1
    1726             : 
    1727   209983032 :          do ig = 1, ng8
    1728   558658752 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1729   558658752 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1730   372439168 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1731   372439168 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1732   372439168 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1733   372439168 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1734             :             abso3 =  (ka_mo3(indminor(lay),ig) + minorfrac(lay) * &
    1735   186219584 :                  (ka_mo3(indminor(lay)+1,ig) - ka_mo3(indminor(lay),ig)))
    1736             :             absn2o =  (ka_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1737   186219584 :                  (ka_mn2o(indminor(lay)+1,ig) - ka_mn2o(indminor(lay),ig)))
    1738   186219584 :             taug(lay,ngs7+ig) = colh2o(lay) * &
    1739   372439168 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    1740   372439168 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1741   372439168 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1742   372439168 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    1743             :                  + tauself + taufor &
    1744             :                  + adjcolco2(lay)*absco2 &
    1745   186219584 :                  + colo3(lay) * abso3 &
    1746   186219584 :                  + coln2o(lay) * absn2o &
    1747   186219584 :                  + wx(3,lay) * cfc12(ig) &
    1748   931097920 :                  + wx(4,lay) * cfc22adj(ig)
    1749   209497032 :             fracs(lay,ngs7+ig) = fracrefa(ig)
    1750             :          enddo
    1751             :       enddo
    1752             : 
    1753             : ! Upper atmosphere loop
    1754    22406552 :       do lay = laytrop+1, nlayers
    1755             : 
    1756             : !  In atmospheres where the amount of CO2 is too great to be considered
    1757             : !  a minor species, adjust the column amount of CO2 by an empirical factor 
    1758             : !  to obtain the proper contribution.
    1759    21920552 :          chi_co2(lay) = colco2(lay)/coldry(lay)
    1760    21920552 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1761    21920552 :          if (ratco2(lay) .gt. 3.0_r8) then
    1762           0 :             adjfac(lay) = 2.0_r8+(ratco2(lay)-2.0_r8)**0.65_r8
    1763           0 :             adjcolco2(lay) = adjfac(lay)*chi_mls(2,jp(lay)+1) * coldry(lay)*1.e-20_r8
    1764             :          else
    1765    21920552 :             adjcolco2(lay) = colco2(lay)
    1766             :          endif
    1767             : 
    1768    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(8) + 1
    1769    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(8) + 1
    1770             : 
    1771   197770968 :          do ig = 1, ng8
    1772   350728832 :             absco2 =  (kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1773   526093248 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig)))
    1774             :             absn2o =  (kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1775   175364416 :                  (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig)))
    1776   175364416 :             taug(lay,ngs7+ig) = colo3(lay) * &
    1777   350728832 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1778   350728832 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1779   350728832 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1780   350728832 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1781             :                  + adjcolco2(lay)*absco2 &
    1782   175364416 :                  + coln2o(lay)*absn2o & 
    1783   175364416 :                  + wx(3,lay) * cfc12(ig) &
    1784   876822080 :                  + wx(4,lay) * cfc22adj(ig)
    1785   197284968 :             fracs(lay,ngs7+ig) = fracrefb(ig)
    1786             :          enddo
    1787             :       enddo
    1788             : 
    1789      486000 :       end subroutine taugb8
    1790             : 
    1791             : !----------------------------------------------------------------------------
    1792      486000 :       subroutine taugb9
    1793             : !----------------------------------------------------------------------------
    1794             : !
    1795             : !     band 9:  1180-1390 cm-1 (low key - h2o,ch4; low minor - n2o)
    1796             : !                             (high key - ch4; high minor - n2o)
    1797             : !----------------------------------------------------------------------------
    1798             : 
    1799             : ! ------- Modules -------
    1800             : 
    1801             :       use parrrtm, only : ng9, ngs8
    1802             :       use rrlw_ref, only : chi_mls
    1803             :       use rrlw_kg09, only : fracrefa, fracrefb, absa, absb, &
    1804             :                             ka_mn2o, kb_mn2o, selfref, forref
    1805             : 
    1806             : ! ------- Declarations -------
    1807             : 
    1808             : ! Local 
    1809      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1810      972000 :       integer, dimension(nlayers) :: js, js1, jmn2o, jpl
    1811      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1812      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1813      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2o, specparm_mn2o, specmult_mn2o, fmn2o
    1814      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1815             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1816      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1817      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    1818             :       real(kind=r8) :: tauself, taufor, n2om1, n2om2, absn2o
    1819      972000 :       real(kind=r8), dimension(nlayers) :: chi_n2o, ratn2o, adjfac, adjcoln2o
    1820             :       real(kind=r8) :: refrat_planck_a, refrat_m_a
    1821             :       real(kind=r8) :: tau_major, tau_major1
    1822             : 
    1823             : 
    1824             : ! Minor gas mapping level :
    1825             : !     lower - n2o, p = 706.272 mbar, t = 278.94 k
    1826             : !     upper - n2o, p = 95.58 mbar, t = 215.7 k
    1827             : 
    1828             : ! Calculate reference ratio to be used in calculation of Planck
    1829             : ! fraction in lower/upper atmosphere.
    1830             : 
    1831             : ! P = 212 mb
    1832      486000 :       refrat_planck_a = chi_mls(1,9)/chi_mls(6,9)
    1833             : 
    1834             : ! P = 706.272 mb 
    1835      486000 :       refrat_m_a = chi_mls(1,3)/chi_mls(6,3)
    1836             : 
    1837             : ! Compute the optical depth by interpolating in ln(pressure), 
    1838             : ! temperature, and appropriate species.  Below laytrop, the water
    1839             : ! vapor self-continuum and foreign continuum is interpolated 
    1840             : ! (in temperature) separately.  
    1841             : 
    1842             : ! Lower atmosphere loop
    1843    23763448 :       do lay = 1, laytrop
    1844             : 
    1845    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    1846    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1847    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1848    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    1849    23277448 :          js(lay) = 1 + int(specmult(lay))
    1850    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1851             : 
    1852    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    1853    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1854    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1855    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1856    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    1857    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1858             : 
    1859    23277448 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colch4(lay)
    1860    23277448 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
    1861    23277448 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
    1862    23277448 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
    1863    23277448 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
    1864    23277448 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
    1865             : 
    1866             : !  In atmospheres where the amount of N2O is too great to be considered
    1867             : !  a minor species, adjust the column amount of N2O by an empirical factor 
    1868             : !  to obtain the proper contribution.
    1869    23277448 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    1870    23277448 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    1871    23277448 :          if (ratn2o(lay) .gt. 1.5_r8) then
    1872           0 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
    1873           0 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
    1874             :          else
    1875    23277448 :             adjcoln2o(lay) = coln2o(lay)
    1876             :          endif
    1877             : 
    1878    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    1879    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1880    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1881    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1882    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1883    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1884             : 
    1885    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(9) + js(lay)
    1886    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(9) + js1(lay)
    1887             : 
    1888    23277448 :          if (specparm(lay) .lt. 0.125_r8) then
    1889     3886970 :             p = fs(lay) - 1
    1890     3886970 :             p4 = p**4
    1891     3886970 :             fk0 = p4
    1892     3886970 :             fk1 = 1 - p - 2.0_r8*p4
    1893     3886970 :             fk2 = p + p4
    1894     3886970 :             fac000(lay) = fk0*fac00(lay)
    1895     3886970 :             fac100(lay) = fk1*fac00(lay)
    1896     3886970 :             fac200(lay) = fk2*fac00(lay)
    1897     3886970 :             fac010(lay) = fk0*fac10(lay)
    1898     3886970 :             fac110(lay) = fk1*fac10(lay)
    1899     3886970 :             fac210(lay) = fk2*fac10(lay)
    1900    19390478 :          else if (specparm(lay) .gt. 0.875_r8) then
    1901       14858 :             p = -fs(lay) 
    1902       14858 :             p4 = p**4
    1903       14858 :             fk0 = p4
    1904       14858 :             fk1 = 1 - p - 2.0_r8*p4
    1905       14858 :             fk2 = p + p4
    1906       14858 :             fac000(lay) = fk0*fac00(lay)
    1907       14858 :             fac100(lay) = fk1*fac00(lay)
    1908       14858 :             fac200(lay) = fk2*fac00(lay)
    1909       14858 :             fac010(lay) = fk0*fac10(lay)
    1910       14858 :             fac110(lay) = fk1*fac10(lay)
    1911       14858 :             fac210(lay) = fk2*fac10(lay)
    1912             :          else
    1913    19375620 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1914    19375620 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1915    19375620 :             fac100(lay) = fs(lay) * fac00(lay)
    1916    19375620 :             fac110(lay) = fs(lay) * fac10(lay)
    1917             :          endif
    1918             : 
    1919    23277448 :          if (specparm1(lay) .lt. 0.125_r8) then
    1920     1558856 :             p = fs1(lay) - 1
    1921     1558856 :             p4 = p**4
    1922     1558856 :             fk0 = p4
    1923     1558856 :             fk1 = 1 - p - 2.0_r8*p4
    1924     1558856 :             fk2 = p + p4
    1925     1558856 :             fac001(lay) = fk0*fac01(lay)
    1926     1558856 :             fac101(lay) = fk1*fac01(lay)
    1927     1558856 :             fac201(lay) = fk2*fac01(lay)
    1928     1558856 :             fac011(lay) = fk0*fac11(lay)
    1929     1558856 :             fac111(lay) = fk1*fac11(lay)
    1930     1558856 :             fac211(lay) = fk2*fac11(lay)
    1931    21718592 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1932      514158 :             p = -fs1(lay) 
    1933      514158 :             p4 = p**4
    1934      514158 :             fk0 = p4
    1935      514158 :             fk1 = 1 - p - 2.0_r8*p4
    1936      514158 :             fk2 = p + p4
    1937      514158 :             fac001(lay) = fk0*fac01(lay)
    1938      514158 :             fac101(lay) = fk1*fac01(lay)
    1939      514158 :             fac201(lay) = fk2*fac01(lay)
    1940      514158 :             fac011(lay) = fk0*fac11(lay)
    1941      514158 :             fac111(lay) = fk1*fac11(lay)
    1942      514158 :             fac211(lay) = fk2*fac11(lay)
    1943             :          else
    1944    21204434 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1945    21204434 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1946    21204434 :             fac101(lay) = fs1(lay) * fac01(lay)
    1947    21204434 :             fac111(lay) = fs1(lay) * fac11(lay)
    1948             :          endif
    1949             : 
    1950   303092824 :          do ig = 1, ng9
    1951   837988128 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1952   837988128 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1953   558658752 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1954   558658752 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1955   837988128 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
    1956   837988128 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
    1957   279329376 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
    1958   279329376 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
    1959   279329376 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
    1960             : 
    1961   279329376 :             if (specparm(lay) .lt. 0.125_r8) then
    1962             :                tau_major = speccomb(lay) * &
    1963    46643640 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1964    46643640 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1965    46643640 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1966    46643640 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1967    46643640 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1968   279861840 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1969   232685736 :             else if (specparm(lay) .gt. 0.875_r8) then
    1970             :                tau_major = speccomb(lay) * &
    1971      178296 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1972      178296 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1973      178296 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1974      178296 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1975      178296 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1976     1069776 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1977             :             else
    1978             :                tau_major = speccomb(lay) * &
    1979   232507440 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1980   232507440 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1981   232507440 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1982   930029760 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1983             :             endif
    1984             : 
    1985   279329376 :             if (specparm1(lay) .lt. 0.125_r8) then
    1986             :                tau_major1 = speccomb1(lay) * &
    1987    18706272 :                     (fac001(lay) * absa(ind1(lay),ig) + & 
    1988    18706272 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1989    18706272 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1990    18706272 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1991    18706272 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1992   112237632 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1993   260623104 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1994             :                tau_major1 = speccomb1(lay) * &
    1995     6169896 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1996     6169896 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1997     6169896 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1998     6169896 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1999     6169896 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2000    37019376 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2001             :             else
    2002             :                tau_major1 = speccomb1(lay) * &
    2003   254453208 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2004   254453208 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2005   254453208 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2006  1017812832 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2007             :             endif
    2008             : 
    2009   279329376 :             taug(lay,ngs8+ig) = tau_major + tau_major1 &
    2010             :                  + tauself + taufor &
    2011   279329376 :                  + adjcoln2o(lay)*absn2o
    2012   558658752 :             fracs(lay,ngs8+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2013   581936200 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2014             :          enddo
    2015             :       enddo
    2016             : 
    2017             : ! Upper atmosphere loop
    2018    22406552 :       do lay = laytrop+1, nlayers
    2019             : 
    2020             : !  In atmospheres where the amount of N2O is too great to be considered
    2021             : !  a minor species, adjust the column amount of N2O by an empirical factor 
    2022             : !  to obtain the proper contribution.
    2023    21920552 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    2024    21920552 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    2025    21920552 :          if (ratn2o(lay) .gt. 1.5_r8) then
    2026    12971441 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
    2027    12971441 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
    2028             :          else
    2029     8949111 :             adjcoln2o(lay) = coln2o(lay)
    2030             :          endif
    2031             : 
    2032    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(9) + 1
    2033    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(9) + 1
    2034             : 
    2035   285453176 :          do ig = 1, ng9
    2036   526093248 :             absn2o = kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    2037   789139872 :                 (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig))
    2038   263046624 :             taug(lay,ngs8+ig) = colch4(lay) * &
    2039   526093248 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2040   526093248 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2041   526093248 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2042   526093248 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2043  1315233120 :                  + adjcoln2o(lay)*absn2o
    2044   284967176 :             fracs(lay,ngs8+ig) = fracrefb(ig)
    2045             :          enddo
    2046             :       enddo
    2047             : 
    2048      486000 :       end subroutine taugb9
    2049             : 
    2050             : !----------------------------------------------------------------------------
    2051      486000 :       subroutine taugb10
    2052             : !----------------------------------------------------------------------------
    2053             : !
    2054             : !     band 10:  1390-1480 cm-1 (low key - h2o; high key - h2o)
    2055             : !----------------------------------------------------------------------------
    2056             : 
    2057             : ! ------- Modules -------
    2058             : 
    2059             :       use parrrtm, only : ng10, ngs9
    2060             :       use rrlw_kg10, only : fracrefa, fracrefb, absa, absb, &
    2061             :                             selfref, forref
    2062             : 
    2063             : ! ------- Declarations -------
    2064             : 
    2065             : ! Local 
    2066      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2067             :       real(kind=r8) :: tauself, taufor
    2068             : 
    2069             : 
    2070             : ! Compute the optical depth by interpolating in ln(pressure) and 
    2071             : ! temperature.  Below laytrop, the water vapor self-continuum and
    2072             : ! foreign continuum is interpolated (in temperature) separately.
    2073             : 
    2074             : ! Lower atmosphere loop
    2075    23763448 :       do lay = 1, laytrop
    2076    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(10) + 1
    2077    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(10) + 1
    2078             : 
    2079   163428136 :          do ig = 1, ng10
    2080   418994064 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2081   418994064 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2082   279329376 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2083   279329376 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2084   139664688 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2085   279329376 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2086   279329376 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2087   279329376 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2088   279329376 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    2089   698323440 :                  + tauself + taufor
    2090   162942136 :             fracs(lay,ngs9+ig) = fracrefa(ig)
    2091             :          enddo
    2092             :       enddo
    2093             : 
    2094             : ! Upper atmosphere loop
    2095    22406552 :       do lay = laytrop+1, nlayers
    2096    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(10) + 1
    2097    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(10) + 1
    2098             : 
    2099   153929864 :          do ig = 1, ng10
    2100   394569936 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2101   394569936 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2102   131523312 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2103   263046624 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2104   263046624 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2105   263046624 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2106   263046624 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2107   657616560 :                  + taufor
    2108   153443864 :             fracs(lay,ngs9+ig) = fracrefb(ig)
    2109             :          enddo
    2110             :       enddo
    2111             : 
    2112      486000 :       end subroutine taugb10
    2113             : 
    2114             : !----------------------------------------------------------------------------
    2115      486000 :       subroutine taugb11
    2116             : !----------------------------------------------------------------------------
    2117             : !
    2118             : !     band 11:  1480-1800 cm-1 (low - h2o; low minor - o2)
    2119             : !                              (high key - h2o; high minor - o2)
    2120             : !----------------------------------------------------------------------------
    2121             : 
    2122             : ! ------- Modules -------
    2123             : 
    2124             :       use parrrtm, only : ng11, ngs10
    2125             :       use rrlw_kg11, only : fracrefa, fracrefb, absa, absb, &
    2126             :                             ka_mo2, kb_mo2, selfref, forref
    2127             : 
    2128             : ! ------- Declarations -------
    2129             : 
    2130             : ! Local 
    2131      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2132      972000 :       real(kind=r8) :: scaleo2(nlayers), tauself, taufor, tauo2
    2133             : 
    2134             : 
    2135             : ! Minor gas mapping level :
    2136             : !     lower - o2, p = 706.2720 mbar, t = 278.94 k
    2137             : !     upper - o2, p = 4.758820 mbarm t = 250.85 k
    2138             : 
    2139             : ! Compute the optical depth by interpolating in ln(pressure) and 
    2140             : ! temperature.  Below laytrop, the water vapor self-continuum and
    2141             : ! foreign continuum is interpolated (in temperature) separately.
    2142             : 
    2143             : ! Lower atmosphere loop
    2144    23763448 :       do lay = 1, laytrop
    2145    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(11) + 1
    2146    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(11) + 1
    2147    23277448 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2148   209983032 :          do ig = 1, ng11
    2149   558658752 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2150   558658752 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2151   372439168 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2152   372439168 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    2153   372439168 :             tauo2 =  scaleo2(lay) * (ka_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2154   372439168 :                  (ka_mo2(indminor(lay)+1,ig) - ka_mo2(indminor(lay),ig)))
    2155   186219584 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2156   372439168 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2157   372439168 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2158   372439168 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2159   372439168 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2160             :                  + tauself + taufor &
    2161   931097920 :                  + tauo2
    2162   209497032 :             fracs(lay,ngs10+ig) = fracrefa(ig)
    2163             :          enddo
    2164             :       enddo
    2165             : 
    2166             : ! Upper atmosphere loop
    2167    22406552 :       do lay = laytrop+1, nlayers
    2168    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(11) + 1
    2169    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(11) + 1
    2170    21920552 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2171   197770968 :          do ig = 1, ng11
    2172   526093248 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2173   526093248 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2174   350728832 :             tauo2 =  scaleo2(lay) * (kb_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2175   350728832 :                  (kb_mo2(indminor(lay)+1,ig) - kb_mo2(indminor(lay),ig)))
    2176   175364416 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2177   350728832 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2178   350728832 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2179   350728832 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2180   350728832 :                  fac11(lay) * absb(ind1(lay)+1,ig))  &
    2181             :                  + taufor &
    2182   876822080 :                  + tauo2
    2183   197284968 :             fracs(lay,ngs10+ig) = fracrefb(ig)
    2184             :          enddo
    2185             :       enddo
    2186             : 
    2187      486000 :       end subroutine taugb11
    2188             : 
    2189             : !----------------------------------------------------------------------------
    2190      486000 :       subroutine taugb12
    2191             : !----------------------------------------------------------------------------
    2192             : !
    2193             : !     band 12:  1800-2080 cm-1 (low - h2o,co2; high - nothing)
    2194             : !----------------------------------------------------------------------------
    2195             : 
    2196             : ! ------- Modules -------
    2197             : 
    2198             :       use parrrtm, only : ng12, ngs11
    2199             :       use rrlw_ref, only : chi_mls
    2200             :       use rrlw_kg12, only : fracrefa, absa, &
    2201             :                             selfref, forref
    2202             : 
    2203             : ! ------- Declarations -------
    2204             : 
    2205             : ! Local 
    2206      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2207      972000 :       integer, dimension(nlayers) :: js, js1, jpl
    2208      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2209      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2210      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2211             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2212      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2213      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    2214             :       real(kind=r8) :: tauself, taufor
    2215             :       real(kind=r8) :: refrat_planck_a
    2216             :       real(kind=r8) :: tau_major, tau_major1
    2217             : 
    2218             : 
    2219             : ! Calculate reference ratio to be used in calculation of Planck
    2220             : ! fraction in lower/upper atmosphere.
    2221             : 
    2222             : ! P =   174.164 mb 
    2223      486000 :       refrat_planck_a = chi_mls(1,10)/chi_mls(2,10)
    2224             : 
    2225             : ! Compute the optical depth by interpolating in ln(pressure), 
    2226             : ! temperature, and appropriate species.  Below laytrop, the water
    2227             : ! vapor self-continuum adn foreign continuum is interpolated 
    2228             : ! (in temperature) separately.  
    2229             : 
    2230             : ! Lower atmosphere loop
    2231    23763448 :       do lay = 1, laytrop
    2232             : 
    2233    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    2234    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2235    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2236    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    2237    23277448 :          js(lay) = 1 + int(specmult(lay))
    2238    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2239             : 
    2240    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    2241    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2242    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2243    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2244    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    2245    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2246             : 
    2247    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    2248    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2249    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2250    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2251    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2252    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2253             : 
    2254    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(12) + js(lay)
    2255    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(12) + js1(lay)
    2256             :       enddo
    2257    23763448 :       do lay = 1, laytrop
    2258    23763448 :          if (specparm(lay) .lt. 0.125_r8) then
    2259     3737204 :             p = fs(lay) - 1
    2260     3737204 :             p4 = p**4
    2261     3737204 :             fk0 = p4
    2262     3737204 :             fk1 = 1 - p - 2.0_r8*p4
    2263     3737204 :             fk2 = p + p4
    2264     3737204 :             fac000(lay) = fk0*fac00(lay)
    2265     3737204 :             fac100(lay) = fk1*fac00(lay)
    2266     3737204 :             fac200(lay) = fk2*fac00(lay)
    2267     3737204 :             fac010(lay) = fk0*fac10(lay)
    2268     3737204 :             fac110(lay) = fk1*fac10(lay)
    2269     3737204 :             fac210(lay) = fk2*fac10(lay)
    2270    19540244 :          else if (specparm(lay) .gt. 0.875_r8) then
    2271       31557 :             p = -fs(lay) 
    2272       31557 :             p4 = p**4
    2273       31557 :             fk0 = p4
    2274       31557 :             fk1 = 1 - p - 2.0_r8*p4
    2275       31557 :             fk2 = p + p4
    2276       31557 :             fac000(lay) = fk0*fac00(lay)
    2277       31557 :             fac100(lay) = fk1*fac00(lay)
    2278       31557 :             fac200(lay) = fk2*fac00(lay)
    2279       31557 :             fac010(lay) = fk0*fac10(lay)
    2280       31557 :             fac110(lay) = fk1*fac10(lay)
    2281       31557 :             fac210(lay) = fk2*fac10(lay)
    2282             :          else
    2283    19508687 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2284    19508687 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2285    19508687 :             fac100(lay) = fs(lay) * fac00(lay)
    2286    19508687 :             fac110(lay) = fs(lay) * fac10(lay)
    2287             :          endif
    2288             : 
    2289             :       enddo
    2290    23763448 :       do lay = 1, laytrop
    2291    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
    2292     1487462 :             p = fs1(lay) - 1
    2293     1487462 :             p4 = p**4
    2294     1487462 :             fk0 = p4
    2295     1487462 :             fk1 = 1 - p - 2.0_r8*p4
    2296     1487462 :             fk2 = p + p4
    2297     1487462 :             fac001(lay) = fk0*fac01(lay)
    2298     1487462 :             fac101(lay) = fk1*fac01(lay)
    2299     1487462 :             fac201(lay) = fk2*fac01(lay)
    2300     1487462 :             fac011(lay) = fk0*fac11(lay)
    2301     1487462 :             fac111(lay) = fk1*fac11(lay)
    2302     1487462 :             fac211(lay) = fk2*fac11(lay)
    2303    21789986 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2304      675028 :             p = -fs1(lay) 
    2305      675028 :             p4 = p**4
    2306      675028 :             fk0 = p4
    2307      675028 :             fk1 = 1 - p - 2.0_r8*p4
    2308      675028 :             fk2 = p + p4
    2309      675028 :             fac001(lay) = fk0*fac01(lay)
    2310      675028 :             fac101(lay) = fk1*fac01(lay)
    2311      675028 :             fac201(lay) = fk2*fac01(lay)
    2312      675028 :             fac011(lay) = fk0*fac11(lay)
    2313      675028 :             fac111(lay) = fk1*fac11(lay)
    2314      675028 :             fac211(lay) = fk2*fac11(lay)
    2315             :          else
    2316    21114958 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2317    21114958 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2318    21114958 :             fac101(lay) = fs1(lay) * fac01(lay)
    2319    21114958 :             fac111(lay) = fs1(lay) * fac11(lay)
    2320             :          endif
    2321             :       enddo
    2322    23763448 :       do lay = 1, laytrop
    2323             : 
    2324   209983032 :          do ig = 1, ng12
    2325   558658752 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2326   744878336 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2327   372439168 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2328   372439168 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2329             : 
    2330   186219584 :             if (specparm(lay) .lt. 0.125_r8) then
    2331             :                tau_major = speccomb(lay) * &
    2332    29897632 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2333    29897632 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2334    29897632 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2335    29897632 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2336    29897632 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2337   179385792 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2338   156321952 :             else if (specparm(lay) .gt. 0.875_r8) then
    2339             :                tau_major = speccomb(lay) * &
    2340      252456 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2341      252456 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2342      252456 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2343      252456 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2344      252456 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2345     1514736 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2346             :             else
    2347             :                tau_major = speccomb(lay) * &
    2348   156069496 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2349   156069496 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2350   156069496 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2351   624277984 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2352             :             endif
    2353             : 
    2354   186219584 :             if (specparm1(lay) .lt. 0.125_r8) then
    2355             :                tau_major1 = speccomb1(lay) * &
    2356    11899696 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2357    11899696 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2358    11899696 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2359    11899696 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2360    11899696 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2361    71398176 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2362   174319888 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2363             :                tau_major1 = speccomb1(lay) * &
    2364     5400224 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2365     5400224 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2366     5400224 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2367     5400224 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2368     5400224 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2369    32401344 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2370             :             else
    2371             :                tau_major1 = speccomb1(lay) * &
    2372   168919664 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2373   168919664 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2374   168919664 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2375   675678656 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2376             :             endif
    2377             : 
    2378   186219584 :             taug(lay,ngs11+ig) = tau_major + tau_major1 &
    2379   186219584 :                  + tauself + taufor
    2380   372439168 :             fracs(lay,ngs11+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2381   395716616 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2382             :          enddo
    2383             :       enddo
    2384             : 
    2385             : ! Upper atmosphere loop
    2386    22406552 :       do lay = laytrop+1, nlayers
    2387   197770968 :          do ig = 1, ng12
    2388   175364416 :             taug(lay,ngs11+ig) = 0.0_r8
    2389   197284968 :             fracs(lay,ngs11+ig) = 0.0_r8
    2390             :          enddo
    2391             :       enddo
    2392             : 
    2393      486000 :       end subroutine taugb12
    2394             : 
    2395             : !----------------------------------------------------------------------------
    2396      486000 :       subroutine taugb13
    2397             : !----------------------------------------------------------------------------
    2398             : !
    2399             : !     band 13:  2080-2250 cm-1 (low key - h2o,n2o; high minor - o3 minor)
    2400             : !----------------------------------------------------------------------------
    2401             : 
    2402             : ! ------- Modules -------
    2403             : 
    2404             :       use parrrtm, only : ng13, ngs12
    2405             :       use rrlw_ref, only : chi_mls
    2406             :       use rrlw_kg13, only : fracrefa, fracrefb, absa, &
    2407             :                             ka_mco2, ka_mco, kb_mo3, selfref, forref
    2408             : 
    2409             : ! ------- Declarations -------
    2410             : 
    2411             : ! Local 
    2412      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2413      972000 :       integer, dimension(nlayers) :: js, js1, jmco2, jmco, jpl
    2414      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2415      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2416      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mco2, specparm_mco2, specmult_mco2, fmco2
    2417      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mco, specparm_mco, specmult_mco, fmco
    2418      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2419             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2420      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2421      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    2422             :       real(kind=r8) :: tauself, taufor, co2m1, co2m2, absco2 
    2423             :       real(kind=r8) :: com1, com2, absco, abso3
    2424      972000 :       real(kind=r8), dimension(nlayers) :: chi_co2, ratco2, adjfac, adjcolco2
    2425             :       real(kind=r8) :: refrat_planck_a, refrat_m_a, refrat_m_a3
    2426             :       real(kind=r8) :: tau_major, tau_major1
    2427             : 
    2428             : 
    2429             : ! Minor gas mapping levels :
    2430             : !     lower - co2, p = 1053.63 mb, t = 294.2 k
    2431             : !     lower - co, p = 706 mb, t = 278.94 k
    2432             : !     upper - o3, p = 95.5835 mb, t = 215.7 k
    2433             : 
    2434             : ! Calculate reference ratio to be used in calculation of Planck
    2435             : ! fraction in lower/upper atmosphere.
    2436             : 
    2437             : ! P = 473.420 mb (Level 5)
    2438      486000 :       refrat_planck_a = chi_mls(1,5)/chi_mls(4,5)
    2439             : 
    2440             : ! P = 1053. (Level 1)
    2441      486000 :       refrat_m_a = chi_mls(1,1)/chi_mls(4,1)
    2442             : 
    2443             : ! P = 706. (Level 3)
    2444      486000 :       refrat_m_a3 = chi_mls(1,3)/chi_mls(4,3)
    2445             : 
    2446             : ! Compute the optical depth by interpolating in ln(pressure), 
    2447             : ! temperature, and appropriate species.  Below laytrop, the water
    2448             : ! vapor self-continuum and foreign continuum is interpolated 
    2449             : ! (in temperature) separately.  
    2450             : 
    2451             : ! Lower atmosphere loop
    2452    23763448 :       do lay = 1, laytrop
    2453             : 
    2454    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2on2o(lay)*coln2o(lay)
    2455    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2456    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2457    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    2458    23277448 :          js(lay) = 1 + int(specmult(lay))
    2459    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2460             : 
    2461    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2on2o_1(lay)*coln2o(lay)
    2462    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2463    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2464    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2465    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    2466    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2467             : 
    2468    23277448 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*coln2o(lay)
    2469    23277448 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    2470    23277448 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    2471    23277448 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    2472    23277448 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    2473    23277448 :          fmco2(lay) = mod(specmult_mco2(lay),1.0_r8)
    2474             : 
    2475             : !  In atmospheres where the amount of CO2 is too great to be considered
    2476             : !  a minor species, adjust the column amount of CO2 by an empirical factor 
    2477             : !  to obtain the proper contribution.
    2478    23277448 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    2479    23763448 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/3.55e-4_r8
    2480             :       enddo
    2481    23763448 :       do lay = 1, laytrop
    2482    23763448 :          if (ratco2(lay) .gt. 3.0_r8) then
    2483           0 :             adjfac(lay) = 2.0_r8+(ratco2(lay)-2.0_r8)**0.68_r8
    2484           0 :             adjcolco2(lay) = adjfac(lay)*3.55e-4*coldry(lay)*1.e-20_r8
    2485             :          else
    2486    23277448 :             adjcolco2(lay) = colco2(lay)
    2487             :          endif
    2488             : 
    2489             :       enddo
    2490    23763448 :       do lay = 1, laytrop
    2491    23277448 :          speccomb_mco(lay) = colh2o(lay) + refrat_m_a3*coln2o(lay)
    2492    23277448 :          specparm_mco(lay) = colh2o(lay)/speccomb_mco(lay)
    2493    23277448 :          if (specparm_mco(lay) .ge. oneminus) specparm_mco(lay) = oneminus
    2494    23277448 :          specmult_mco(lay) = 8._r8*specparm_mco(lay)
    2495    23277448 :          jmco(lay) = 1 + int(specmult_mco(lay))
    2496    23277448 :          fmco(lay) = mod(specmult_mco(lay),1.0_r8)
    2497             : 
    2498    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*coln2o(lay)
    2499    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2500    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2501    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2502    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2503    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2504             : 
    2505    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(13) + js(lay)
    2506    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(13) + js1(lay)
    2507             :       enddo
    2508    23763448 :       do lay = 1, laytrop
    2509             : 
    2510    23763448 :          if (specparm(lay) .lt. 0.125_r8) then
    2511     3650444 :             p = fs(lay) - 1
    2512     3650444 :             p4 = p**4
    2513     3650444 :             fk0 = p4
    2514     3650444 :             fk1 = 1 - p - 2.0_r8*p4
    2515     3650444 :             fk2 = p + p4
    2516     3650444 :             fac000(lay) = fk0*fac00(lay)
    2517     3650444 :             fac100(lay) = fk1*fac00(lay)
    2518     3650444 :             fac200(lay) = fk2*fac00(lay)
    2519     3650444 :             fac010(lay) = fk0*fac10(lay)
    2520     3650444 :             fac110(lay) = fk1*fac10(lay)
    2521     3650444 :             fac210(lay) = fk2*fac10(lay)
    2522    19627004 :          else if (specparm(lay) .gt. 0.875_r8) then
    2523       21290 :             p = -fs(lay) 
    2524       21290 :             p4 = p**4
    2525       21290 :             fk0 = p4
    2526       21290 :             fk1 = 1 - p - 2.0_r8*p4
    2527       21290 :             fk2 = p + p4
    2528       21290 :             fac000(lay) = fk0*fac00(lay)
    2529       21290 :             fac100(lay) = fk1*fac00(lay)
    2530       21290 :             fac200(lay) = fk2*fac00(lay)
    2531       21290 :             fac010(lay) = fk0*fac10(lay)
    2532       21290 :             fac110(lay) = fk1*fac10(lay)
    2533       21290 :             fac210(lay) = fk2*fac10(lay)
    2534             :          else
    2535    19605714 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2536    19605714 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2537    19605714 :             fac100(lay) = fs(lay) * fac00(lay)
    2538    19605714 :             fac110(lay) = fs(lay) * fac10(lay)
    2539             :          endif
    2540             : 
    2541             :       enddo
    2542    23763448 :       do lay = 1, laytrop
    2543    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
    2544     1431721 :             p = fs1(lay) - 1
    2545     1431721 :             p4 = p**4
    2546     1431721 :             fk0 = p4
    2547     1431721 :             fk1 = 1 - p - 2.0_r8*p4
    2548     1431721 :             fk2 = p + p4
    2549     1431721 :             fac001(lay) = fk0*fac01(lay)
    2550     1431721 :             fac101(lay) = fk1*fac01(lay)
    2551     1431721 :             fac201(lay) = fk2*fac01(lay)
    2552     1431721 :             fac011(lay) = fk0*fac11(lay)
    2553     1431721 :             fac111(lay) = fk1*fac11(lay)
    2554     1431721 :             fac211(lay) = fk2*fac11(lay)
    2555    21845727 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2556      607371 :             p = -fs1(lay) 
    2557      607371 :             p4 = p**4
    2558      607371 :             fk0 = p4
    2559      607371 :             fk1 = 1 - p - 2.0_r8*p4
    2560      607371 :             fk2 = p + p4
    2561      607371 :             fac001(lay) = fk0*fac01(lay)
    2562      607371 :             fac101(lay) = fk1*fac01(lay)
    2563      607371 :             fac201(lay) = fk2*fac01(lay)
    2564      607371 :             fac011(lay) = fk0*fac11(lay)
    2565      607371 :             fac111(lay) = fk1*fac11(lay)
    2566      607371 :             fac211(lay) = fk2*fac11(lay)
    2567             :          else
    2568    21238356 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2569    21238356 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2570    21238356 :             fac101(lay) = fs1(lay) * fac01(lay)
    2571    21238356 :             fac111(lay) = fs1(lay) * fac11(lay)
    2572             :          endif
    2573             :       enddo
    2574    23763448 :       do lay = 1, laytrop
    2575             : 
    2576   116873240 :          do ig = 1, ng13
    2577   279329376 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2578   372439168 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2579   186219584 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2580   186219584 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2581   279329376 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    2582   279329376 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    2583    93109792 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    2584    93109792 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    2585    93109792 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    2586    93109792 :             com1 = ka_mco(jmco(lay),indminor(lay),ig) + fmco(lay) * &
    2587   186219584 :                  (ka_mco(jmco(lay)+1,indminor(lay),ig) - ka_mco(jmco(lay),indminor(lay),ig))
    2588             :             com2 = ka_mco(jmco(lay),indminor(lay)+1,ig) + fmco(lay) * &
    2589    93109792 :                  (ka_mco(jmco(lay)+1,indminor(lay)+1,ig) - ka_mco(jmco(lay),indminor(lay)+1,ig))
    2590    93109792 :             absco = com1 + minorfrac(lay) * (com2 - com1)
    2591             : 
    2592    93109792 :             if (specparm(lay) .lt. 0.125_r8) then
    2593             :                tau_major = speccomb(lay) * &
    2594    14601776 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2595    14601776 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2596    14601776 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2597    14601776 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2598    14601776 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2599    87610656 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2600    78508016 :             else if (specparm(lay) .gt. 0.875_r8) then
    2601             :                tau_major = speccomb(lay) * &
    2602       85160 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2603       85160 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2604       85160 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2605       85160 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2606       85160 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2607      510960 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2608             :             else
    2609             :                tau_major = speccomb(lay) * &
    2610    78422856 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2611    78422856 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2612    78422856 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2613   313691424 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2614             :             endif
    2615             : 
    2616    93109792 :             if (specparm1(lay) .lt. 0.125_r8) then
    2617             :                tau_major1 = speccomb1(lay) * &
    2618     5726884 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2619     5726884 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2620     5726884 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2621     5726884 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2622     5726884 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2623    34361304 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2624    87382908 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2625             :                tau_major1 = speccomb1(lay) * &
    2626     2429484 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2627     2429484 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2628     2429484 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2629     2429484 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2630     2429484 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2631    14576904 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2632             :             else
    2633             :                tau_major1 = speccomb1(lay) * &
    2634    84953424 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2635    84953424 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2636    84953424 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2637   339813696 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2638             :             endif
    2639             : 
    2640    93109792 :             taug(lay,ngs12+ig) = tau_major + tau_major1 &
    2641             :                  + tauself + taufor &
    2642             :                  + adjcolco2(lay)*absco2 &
    2643   186219584 :                  + colco(lay)*absco
    2644   186219584 :             fracs(lay,ngs12+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2645   209497032 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2646             :          enddo
    2647             :       enddo
    2648             : 
    2649             : ! Upper atmosphere loop
    2650    22406552 :       do lay = laytrop+1, nlayers
    2651   110088760 :          do ig = 1, ng13
    2652   263046624 :             abso3 = kb_mo3(indminor(lay),ig) + minorfrac(lay) * &
    2653   350728832 :                  (kb_mo3(indminor(lay)+1,ig) - kb_mo3(indminor(lay),ig))
    2654    87682208 :             taug(lay,ngs12+ig) = colo3(lay)*abso3
    2655   109602760 :             fracs(lay,ngs12+ig) =  fracrefb(ig)
    2656             :          enddo
    2657             :       enddo
    2658             : 
    2659      486000 :       end subroutine taugb13
    2660             : 
    2661             : !----------------------------------------------------------------------------
    2662      486000 :       subroutine taugb14
    2663             : !----------------------------------------------------------------------------
    2664             : !
    2665             : !     band 14:  2250-2380 cm-1 (low - co2; high - co2)
    2666             : !----------------------------------------------------------------------------
    2667             : 
    2668             : ! ------- Modules -------
    2669             : 
    2670             :       use parrrtm, only : ng14, ngs13
    2671             :       use rrlw_kg14, only : fracrefa, fracrefb, absa, absb, &
    2672             :                             selfref, forref
    2673             : 
    2674             : ! ------- Declarations -------
    2675             : 
    2676             : ! Local 
    2677      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2678             :       real(kind=r8) :: tauself, taufor
    2679             : 
    2680             : 
    2681             : ! Compute the optical depth by interpolating in ln(pressure) and 
    2682             : ! temperature.  Below laytrop, the water vapor self-continuum 
    2683             : ! and foreign continuum is interpolated (in temperature) separately.  
    2684             : 
    2685             : ! Lower atmosphere loop
    2686    23763448 :       do lay = 1, laytrop
    2687    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(14) + 1
    2688    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(14) + 1
    2689    70318344 :          do ig = 1, ng14
    2690   139664688 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2691   139664688 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2692    93109792 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2693    93109792 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2694    46554896 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2695    93109792 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2696    93109792 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2697    93109792 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2698    93109792 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2699   232774480 :                  + tauself + taufor
    2700    69832344 :             fracs(lay,ngs13+ig) = fracrefa(ig)
    2701             :          enddo
    2702             :       enddo
    2703             : 
    2704             : ! Upper atmosphere loop
    2705    22406552 :       do lay = laytrop+1, nlayers
    2706    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(14) + 1
    2707    21920552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(14) + 1
    2708    66247656 :          do ig = 1, ng14
    2709    43841104 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2710   131523312 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2711    87682208 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2712    87682208 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2713   263046624 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    2714    65761656 :             fracs(lay,ngs13+ig) = fracrefb(ig)
    2715             :          enddo
    2716             :       enddo
    2717             : 
    2718      486000 :       end subroutine taugb14
    2719             : 
    2720             : !----------------------------------------------------------------------------
    2721      486000 :       subroutine taugb15
    2722             : !----------------------------------------------------------------------------
    2723             : !
    2724             : !     band 15:  2380-2600 cm-1 (low - n2o,co2; low minor - n2)
    2725             : !                              (high - nothing)
    2726             : !----------------------------------------------------------------------------
    2727             : 
    2728             : ! ------- Modules -------
    2729             : 
    2730             :       use parrrtm, only : ng15, ngs14
    2731             :       use rrlw_ref, only : chi_mls
    2732             :       use rrlw_kg15, only : fracrefa, absa, &
    2733             :                             ka_mn2, selfref, forref
    2734             : 
    2735             : ! ------- Declarations -------
    2736             : 
    2737             : ! Local 
    2738      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2739      972000 :       integer, dimension(nlayers) :: js, js1, jmn2, jpl
    2740      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2741      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2742      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2, specparm_mn2, specmult_mn2, fmn2
    2743      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2744             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2745      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2746      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    2747      972000 :       real(kind=r8) :: scalen2(nlayers), tauself, taufor, n2m1, n2m2, taun2 
    2748             :       real(kind=r8) :: refrat_planck_a, refrat_m_a
    2749             :       real(kind=r8) :: tau_major, tau_major1
    2750             : 
    2751             : 
    2752             : ! Minor gas mapping level : 
    2753             : !     Lower - Nitrogen Continuum, P = 1053., T = 294.
    2754             : 
    2755             : ! Calculate reference ratio to be used in calculation of Planck
    2756             : ! fraction in lower atmosphere.
    2757             : ! P = 1053. mb (Level 1)
    2758      486000 :       refrat_planck_a = chi_mls(4,1)/chi_mls(2,1)
    2759             : 
    2760             : ! P = 1053.
    2761      486000 :       refrat_m_a = chi_mls(4,1)/chi_mls(2,1)
    2762             : 
    2763             : ! Compute the optical depth by interpolating in ln(pressure), 
    2764             : ! temperature, and appropriate species.  Below laytrop, the water
    2765             : ! vapor self-continuum and foreign continuum is interpolated 
    2766             : ! (in temperature) separately.  
    2767             : 
    2768             : ! Lower atmosphere loop
    2769    23763448 :       do lay = 1, laytrop
    2770             : 
    2771    23277448 :          speccomb(lay) = coln2o(lay) + rat_n2oco2(lay)*colco2(lay)
    2772    23277448 :          specparm(lay) = coln2o(lay)/speccomb(lay)
    2773    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2774    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    2775    23277448 :          js(lay) = 1 + int(specmult(lay))
    2776    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2777             : 
    2778    23277448 :          speccomb1(lay) = coln2o(lay) + rat_n2oco2_1(lay)*colco2(lay)
    2779    23277448 :          specparm1(lay) = coln2o(lay)/speccomb1(lay)
    2780    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2781    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2782    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    2783    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2784             : 
    2785    23277448 :          speccomb_mn2(lay) = coln2o(lay) + refrat_m_a*colco2(lay)
    2786    23277448 :          specparm_mn2(lay) = coln2o(lay)/speccomb_mn2(lay)
    2787    23277448 :          if (specparm_mn2(lay) .ge. oneminus) specparm_mn2(lay) = oneminus
    2788    23277448 :          specmult_mn2(lay) = 8._r8*specparm_mn2(lay)
    2789    23277448 :          jmn2(lay) = 1 + int(specmult_mn2(lay))
    2790    23277448 :          fmn2(lay) = mod(specmult_mn2(lay),1.0_r8)
    2791             : 
    2792    23277448 :          speccomb_planck(lay) = coln2o(lay)+refrat_planck_a*colco2(lay)
    2793    23277448 :          specparm_planck(lay) = coln2o(lay)/speccomb_planck(lay)
    2794    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2795    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2796    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2797    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2798             : 
    2799    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(15) + js(lay)
    2800    23277448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(15) + js1(lay)
    2801             :          
    2802    23763448 :          scalen2(lay) = colbrd(lay)*scaleminor(lay)
    2803             :       enddo
    2804    23763448 :       do lay = 1, laytrop
    2805    23763448 :          if (specparm(lay) .lt. 0.125_r8) then
    2806           0 :             p = fs(lay) - 1
    2807           0 :             p4 = p**4
    2808           0 :             fk0 = p4
    2809           0 :             fk1 = 1 - p - 2.0_r8*p4
    2810           0 :             fk2 = p + p4
    2811           0 :             fac000(lay) = fk0*fac00(lay)
    2812           0 :             fac100(lay) = fk1*fac00(lay)
    2813           0 :             fac200(lay) = fk2*fac00(lay)
    2814           0 :             fac010(lay) = fk0*fac10(lay)
    2815           0 :             fac110(lay) = fk1*fac10(lay)
    2816           0 :             fac210(lay) = fk2*fac10(lay)
    2817    23277448 :          else if (specparm(lay) .gt. 0.875_r8) then
    2818           0 :             p = -fs(lay) 
    2819           0 :             p4 = p**4
    2820           0 :             fk0 = p4
    2821           0 :             fk1 = 1 - p - 2.0_r8*p4
    2822           0 :             fk2 = p + p4
    2823           0 :             fac000(lay) = fk0*fac00(lay)
    2824           0 :             fac100(lay) = fk1*fac00(lay)
    2825           0 :             fac200(lay) = fk2*fac00(lay)
    2826           0 :             fac010(lay) = fk0*fac10(lay)
    2827           0 :             fac110(lay) = fk1*fac10(lay)
    2828           0 :             fac210(lay) = fk2*fac10(lay)
    2829             :          else
    2830    23277448 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2831    23277448 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2832    23277448 :             fac100(lay) = fs(lay) * fac00(lay)
    2833    23277448 :             fac110(lay) = fs(lay) * fac10(lay)
    2834             :          endif
    2835             :       enddo
    2836    23763448 :       do lay = 1, laytrop
    2837    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
    2838           0 :             p = fs1(lay) - 1
    2839           0 :             p4 = p**4
    2840           0 :             fk0 = p4
    2841           0 :             fk1 = 1 - p - 2.0_r8*p4
    2842           0 :             fk2 = p + p4
    2843           0 :             fac001(lay) = fk0*fac01(lay)
    2844           0 :             fac101(lay) = fk1*fac01(lay)
    2845           0 :             fac201(lay) = fk2*fac01(lay)
    2846           0 :             fac011(lay) = fk0*fac11(lay)
    2847           0 :             fac111(lay) = fk1*fac11(lay)
    2848           0 :             fac211(lay) = fk2*fac11(lay)
    2849    23277448 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2850           0 :             p = -fs1(lay) 
    2851           0 :             p4 = p**4
    2852           0 :             fk0 = p4
    2853           0 :             fk1 = 1 - p - 2.0_r8*p4
    2854           0 :             fk2 = p + p4
    2855           0 :             fac001(lay) = fk0*fac01(lay)
    2856           0 :             fac101(lay) = fk1*fac01(lay)
    2857           0 :             fac201(lay) = fk2*fac01(lay)
    2858           0 :             fac011(lay) = fk0*fac11(lay)
    2859           0 :             fac111(lay) = fk1*fac11(lay)
    2860           0 :             fac211(lay) = fk2*fac11(lay)
    2861             :          else
    2862    23277448 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2863    23277448 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2864    23277448 :             fac101(lay) = fs1(lay) * fac01(lay)
    2865    23277448 :             fac111(lay) = fs1(lay) * fac11(lay)
    2866             :          endif
    2867             : 
    2868             :       enddo
    2869    23763448 :       do lay = 1, laytrop
    2870    70318344 :          do ig = 1, ng15
    2871   139664688 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2872   186219584 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2873    93109792 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2874    93109792 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2875   139664688 :             n2m1 = ka_mn2(jmn2(lay),indminor(lay),ig) + fmn2(lay) * &
    2876   139664688 :                  (ka_mn2(jmn2(lay)+1,indminor(lay),ig) - ka_mn2(jmn2(lay),indminor(lay),ig))
    2877    46554896 :             n2m2 = ka_mn2(jmn2(lay),indminor(lay)+1,ig) + fmn2(lay) * &
    2878    46554896 :                  (ka_mn2(jmn2(lay)+1,indminor(lay)+1,ig) - ka_mn2(jmn2(lay),indminor(lay)+1,ig))
    2879    46554896 :             taun2 = scalen2(lay) * (n2m1 + minorfrac(lay) * (n2m2 - n2m1))
    2880             : 
    2881    46554896 :             if (specparm(lay) .lt. 0.125_r8) then
    2882             :                tau_major = speccomb(lay) * &
    2883           0 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2884           0 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2885           0 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2886           0 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2887           0 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2888           0 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2889    46554896 :             else if (specparm(lay) .gt. 0.875_r8) then
    2890             :                tau_major = speccomb(lay) * &
    2891           0 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2892           0 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2893           0 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2894           0 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2895           0 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2896           0 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2897             :             else
    2898             :                tau_major = speccomb(lay) * &
    2899    46554896 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2900    46554896 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2901    46554896 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2902   186219584 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2903             :             endif 
    2904             : 
    2905    46554896 :             if (specparm1(lay) .lt. 0.125_r8) then
    2906             :                tau_major1 = speccomb1(lay) * &
    2907           0 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2908           0 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2909           0 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2910           0 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2911           0 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2912           0 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2913    46554896 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2914             :                tau_major1 = speccomb1(lay) * &
    2915           0 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2916           0 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2917           0 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2918           0 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2919           0 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2920           0 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2921             :             else
    2922             :                tau_major1 = speccomb1(lay) * &
    2923    46554896 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2924    46554896 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2925    46554896 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2926   186219584 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2927             :             endif
    2928             : 
    2929    46554896 :             taug(lay,ngs14+ig) = tau_major + tau_major1 &
    2930             :                  + tauself + taufor &
    2931    46554896 :                  + taun2
    2932    93109792 :             fracs(lay,ngs14+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2933   116387240 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2934             :          enddo
    2935             :       enddo
    2936             : 
    2937             : ! Upper atmosphere loop
    2938    22406552 :       do lay = laytrop+1, nlayers
    2939    66247656 :          do ig = 1, ng15
    2940    43841104 :             taug(lay,ngs14+ig) = 0.0_r8
    2941    65761656 :             fracs(lay,ngs14+ig) = 0.0_r8
    2942             :          enddo
    2943             :       enddo
    2944             : 
    2945      486000 :       end subroutine taugb15
    2946             : 
    2947             : !----------------------------------------------------------------------------
    2948      486000 :       subroutine taugb16
    2949             : !----------------------------------------------------------------------------
    2950             : !
    2951             : !     band 16:  2600-3250 cm-1 (low key- h2o,ch4; high key - ch4)
    2952             : !----------------------------------------------------------------------------
    2953             : 
    2954             : ! ------- Modules -------
    2955             : 
    2956             :       use parrrtm, only : ng16, ngs15
    2957             :       use rrlw_ref, only : chi_mls
    2958             :       use rrlw_kg16, only : fracrefa, fracrefb, absa, absb, &
    2959             :                             selfref, forref
    2960             : 
    2961             : ! ------- Declarations -------
    2962             : 
    2963             : ! Local 
    2964      972000 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2965      972000 :       integer, dimension(nlayers) :: js, js1, jpl
    2966      972000 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2967      972000 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2968      972000 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2969             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2970      972000 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2971      972000 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    2972             :       real(kind=r8) :: tauself, taufor
    2973             :       real(kind=r8) :: refrat_planck_a
    2974             :       real(kind=r8) :: tau_major, tau_major1
    2975             : 
    2976             : 
    2977             : ! Calculate reference ratio to be used in calculation of Planck
    2978             : ! fraction in lower atmosphere.
    2979             : 
    2980             : ! P = 387. mb (Level 6)
    2981      486000 :       refrat_planck_a = chi_mls(1,6)/chi_mls(6,6)
    2982             : 
    2983             : ! Compute the optical depth by interpolating in ln(pressure), 
    2984             : ! temperature,and appropriate species.  Below laytrop, the water
    2985             : ! vapor self-continuum and foreign continuum is interpolated 
    2986             : ! (in temperature) separately.  
    2987             : 
    2988             : ! Lower atmosphere loop
    2989    23763448 :       do lay = 1, laytrop
    2990             : 
    2991    23277448 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    2992    23277448 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2993    23277448 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2994    23277448 :          specmult(lay) = 8._r8*(specparm(lay))
    2995    23277448 :          js(lay) = 1 + int(specmult(lay))
    2996    23277448 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2997             : 
    2998    23277448 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    2999    23277448 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    3000    23277448 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    3001    23277448 :          specmult1(lay) = 8._r8*(specparm1(lay))
    3002    23277448 :          js1(lay) = 1 + int(specmult1(lay))
    3003    23277448 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    3004             : 
    3005    23277448 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    3006    23277448 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    3007    23277448 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    3008    23277448 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    3009    23277448 :          jpl(lay)= 1 + int(specmult_planck(lay))
    3010    23277448 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    3011             : 
    3012    23277448 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(16) + js(lay)
    3013    23763448 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(16) + js1(lay)
    3014             :       enddo
    3015    23763448 :       do lay = 1, laytrop
    3016    23763448 :          if (specparm(lay) .lt. 0.125_r8) then
    3017     3886970 :             p = fs(lay) - 1
    3018     3886970 :             p4 = p**4
    3019     3886970 :             fk0 = p4
    3020     3886970 :             fk1 = 1 - p - 2.0_r8*p4
    3021     3886970 :             fk2 = p + p4
    3022     3886970 :             fac000(lay) = fk0*fac00(lay)
    3023     3886970 :             fac100(lay) = fk1*fac00(lay)
    3024     3886970 :             fac200(lay) = fk2*fac00(lay)
    3025     3886970 :             fac010(lay) = fk0*fac10(lay)
    3026     3886970 :             fac110(lay) = fk1*fac10(lay)
    3027     3886970 :             fac210(lay) = fk2*fac10(lay)
    3028    19390478 :          else if (specparm(lay) .gt. 0.875_r8) then
    3029       14858 :             p = -fs(lay) 
    3030       14858 :             p4 = p**4
    3031       14858 :             fk0 = p4
    3032       14858 :             fk1 = 1 - p - 2.0_r8*p4
    3033       14858 :             fk2 = p + p4
    3034       14858 :             fac000(lay) = fk0*fac00(lay)
    3035       14858 :             fac100(lay) = fk1*fac00(lay)
    3036       14858 :             fac200(lay) = fk2*fac00(lay)
    3037       14858 :             fac010(lay) = fk0*fac10(lay)
    3038       14858 :             fac110(lay) = fk1*fac10(lay)
    3039       14858 :             fac210(lay) = fk2*fac10(lay)
    3040             :          else
    3041    19375620 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    3042    19375620 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    3043    19375620 :             fac100(lay) = fs(lay) * fac00(lay)
    3044    19375620 :             fac110(lay) = fs(lay) * fac10(lay)
    3045             :          endif
    3046             :       enddo
    3047    23763448 :       do lay = 1, laytrop
    3048             : 
    3049    23763448 :          if (specparm1(lay) .lt. 0.125_r8) then
    3050     1558856 :             p = fs1(lay) - 1
    3051     1558856 :             p4 = p**4
    3052     1558856 :             fk0 = p4
    3053     1558856 :             fk1 = 1 - p - 2.0_r8*p4
    3054     1558856 :             fk2 = p + p4
    3055     1558856 :             fac001(lay) = fk0*fac01(lay)
    3056     1558856 :             fac101(lay) = fk1*fac01(lay)
    3057     1558856 :             fac201(lay) = fk2*fac01(lay)
    3058     1558856 :             fac011(lay) = fk0*fac11(lay)
    3059     1558856 :             fac111(lay) = fk1*fac11(lay)
    3060     1558856 :             fac211(lay) = fk2*fac11(lay)
    3061    21718592 :          else if (specparm1(lay) .gt. 0.875_r8) then
    3062      514158 :             p = -fs1(lay) 
    3063      514158 :             p4 = p**4
    3064      514158 :             fk0 = p4
    3065      514158 :             fk1 = 1 - p - 2.0_r8*p4
    3066      514158 :             fk2 = p + p4
    3067      514158 :             fac001(lay) = fk0*fac01(lay)
    3068      514158 :             fac101(lay) = fk1*fac01(lay)
    3069      514158 :             fac201(lay) = fk2*fac01(lay)
    3070      514158 :             fac011(lay) = fk0*fac11(lay)
    3071      514158 :             fac111(lay) = fk1*fac11(lay)
    3072      514158 :             fac211(lay) = fk2*fac11(lay)
    3073             :          else
    3074    21204434 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    3075    21204434 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    3076    21204434 :             fac101(lay) = fs1(lay) * fac01(lay)
    3077    21204434 :             fac111(lay) = fs1(lay) * fac11(lay)
    3078             :          endif
    3079             :       enddo
    3080    23763448 :       do lay = 1, laytrop
    3081             : 
    3082    70318344 :          do ig = 1, ng16
    3083   139664688 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    3084   186219584 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    3085    93109792 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    3086    93109792 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    3087             : 
    3088    46554896 :             if (specparm(lay) .lt. 0.125_r8) then
    3089             :                tau_major = speccomb(lay) * &
    3090     7773940 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3091     7773940 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3092     7773940 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    3093     7773940 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3094     7773940 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    3095    46643640 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    3096    38780956 :             else if (specparm(lay) .gt. 0.875_r8) then
    3097             :                tau_major = speccomb(lay) * &
    3098       29716 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    3099       29716 :                     fac100(lay) * absa(ind0(lay),ig) + &
    3100       29716 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    3101       29716 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    3102       29716 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    3103      178296 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    3104             :             else
    3105             :                tau_major = speccomb(lay) * &
    3106    38751240 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3107    38751240 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3108    38751240 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3109   155004960 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    3110             :             endif
    3111             : 
    3112    46554896 :             if (specparm1(lay) .lt. 0.125_r8) then
    3113             :                tau_major1 = speccomb1(lay) * &
    3114     3117712 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3115     3117712 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3116     3117712 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    3117     3117712 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3118     3117712 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    3119    18706272 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    3120    43437184 :             else if (specparm1(lay) .gt. 0.875_r8) then
    3121             :                tau_major1 = speccomb1(lay) * &
    3122     1028316 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    3123     1028316 :                     fac101(lay) * absa(ind1(lay),ig) + &
    3124     1028316 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    3125     1028316 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    3126     1028316 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    3127     6169896 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    3128             :             else
    3129             :                tau_major1 = speccomb1(lay) * &
    3130    42408868 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3131    42408868 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3132    42408868 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3133   169635472 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    3134             :             endif
    3135             : 
    3136    46554896 :             taug(lay,ngs15+ig) = tau_major + tau_major1 &
    3137    46554896 :                  + tauself + taufor
    3138    93109792 :             fracs(lay,ngs15+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    3139   116387240 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    3140             :          enddo
    3141             :       enddo
    3142             : 
    3143             : ! Upper atmosphere loop
    3144    22406552 :       do lay = laytrop+1, nlayers
    3145    21920552 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(16) + 1
    3146    22406552 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(16) + 1
    3147             :       enddo
    3148    22406552 :       do lay = laytrop+1, nlayers
    3149    66247656 :          do ig = 1, ng16
    3150    87682208 :             taug(lay,ngs15+ig) = colch4(lay) * &
    3151   131523312 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    3152    87682208 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    3153    87682208 :                  fac01(lay) * absb(ind1(lay),ig) + &
    3154   306887728 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    3155    65761656 :             fracs(lay,ngs15+ig) = fracrefb(ig)
    3156             :          enddo
    3157             :       enddo
    3158             : 
    3159      486000 :       end subroutine taugb16
    3160             : 
    3161             :       end subroutine taumol
    3162             : 
    3163             :       end module rrtmg_lw_taumol
    3164             : 

Generated by: LCOV version 1.14