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-14 01:18:36 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      552960 :       subroutine taumol(nlayers, pavel, wx, coldry, &
      32      552960 :                         laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
      33      552960 :                         colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
      34      552960 :                         colbrd, fac00, fac01, fac10, fac11, &
      35      552960 :                         rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
      36      552960 :                         rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
      37      552960 :                         rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
      38      552960 :                         selffac, selffrac, indself, forfac, forfrac, indfor, &
      39      552960 :                         minorfrac, scaleminor, scaleminorn2, indminor, &
      40      552960 :                         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      552960 :       call taugb1
     259      552960 :       call taugb2
     260      552960 :       call taugb3
     261      552960 :       call taugb4
     262      552960 :       call taugb5
     263      552960 :       call taugb6
     264      552960 :       call taugb7
     265      552960 :       call taugb8
     266      552960 :       call taugb9
     267      552960 :       call taugb10
     268      552960 :       call taugb11
     269      552960 :       call taugb12
     270      552960 :       call taugb13
     271      552960 :       call taugb14
     272      552960 :       call taugb15
     273      552960 :       call taugb16
     274             : 
     275             :       contains
     276             : 
     277             : !----------------------------------------------------------------------------
     278      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     302     1105920 :       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    19906560 :       do lay = 1, laytrop
     315             : 
     316    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(1) + 1
     317    19353600 :          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    19353600 :          corradj(lay) =  1.
     323    19906560 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     324             :       enddo
     325    19906560 :       do lay = 1, laytrop
     326    19906560 :          if (pavel(lay) .lt. 250._r8) then
     327     3324504 :             corradj(lay) = 1._r8 - 0.15_r8 * (250._r8-pavel(lay)) / 154.4_r8
     328             :          endif
     329             :       enddo
     330    19906560 :       do lay = 1, laytrop
     331   213442560 :          do ig = 1, ng1
     332   580608000 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     333   774144000 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     334   387072000 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     335   387072000 :                  (forref(indfor(lay)+1,ig) -  forref(indfor(lay),ig))) 
     336   387072000 :             taun2 = scalen2(lay)*(ka_mn2(indminor(lay),ig) + & 
     337   387072000 :                  minorfrac(lay) * (ka_mn2(indminor(lay)+1,ig) - ka_mn2(indminor(lay),ig)))
     338   193536000 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     339   387072000 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     340   387072000 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     341   387072000 :                  fac01(lay) * absa(ind1(lay),ig) + &
     342   387072000 :                  fac11(lay) * absa(ind1(lay)+1,ig)) & 
     343   774144000 :                  + tauself + taufor + taun2)
     344   212889600 :              fracs(lay,ig) = fracrefa(ig)
     345             :          enddo
     346             :       enddo
     347             : 
     348             : ! Upper atmosphere loop
     349    12718080 :       do lay = laytrop+1, nlayers
     350             : 
     351    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(1) + 1
     352    12165120 :          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    12165120 :          corradj(lay) =  1._r8 - 0.15_r8 * (pavel(lay) / 95.6_r8)
     357             : 
     358    12718080 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     359             :       enddo
     360    12718080 :       do lay = laytrop+1, nlayers
     361   134369280 :          do ig = 1, ng1
     362   364953600 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     363   486604800 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     364   243302400 :             taun2 = scalen2(lay)*(kb_mn2(indminor(lay),ig) + & 
     365   243302400 :                  minorfrac(lay) * (kb_mn2(indminor(lay)+1,ig) - kb_mn2(indminor(lay),ig)))
     366   121651200 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     367   243302400 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     368   243302400 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     369   243302400 :                  fac01(lay) * absb(ind1(lay),ig) + &
     370   243302400 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &  
     371   486604800 :                  + taufor + taun2)
     372   133816320 :             fracs(lay,ig) = fracrefb(ig)
     373             :          enddo
     374             :       enddo
     375             : 
     376      552960 :       end subroutine taugb1
     377             : 
     378             : !----------------------------------------------------------------------------
     379      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     398     1105920 :       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    19906560 :       do lay = 1, laytrop
     407             : 
     408    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(2) + 1
     409    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(2) + 1
     410             :          !inds = indself(lay)
     411             :          !indf = indfor(lay)
     412             :          !pp = pavel(lay)
     413    19906560 :          corradj(lay) = 1._r8 - .05_r8 * (pavel(lay) - 100._r8) / 900._r8
     414             :       enddo
     415    19906560 :       do lay = 1, laytrop
     416   252149760 :          do ig = 1, ng2
     417   696729600 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     418   928972800 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     419   464486400 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     420   464486400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     421   232243200 :             taug(lay,ngs1+ig) = corradj(lay) * (colh2o(lay) * &
     422   464486400 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     423   464486400 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     424   464486400 :                  fac01(lay) * absa(ind1(lay),ig) + &
     425   464486400 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
     426  1161216000 :                  + tauself + taufor)
     427   251596800 :             fracs(lay,ngs1+ig) = fracrefa(ig)
     428             :          enddo
     429             :       enddo
     430             : 
     431             : ! Upper atmosphere loop
     432    12718080 :       do lay = laytrop+1, nlayers
     433             : 
     434    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(2) + 1
     435    12718080 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(2) + 1
     436             :          !indf = indfor(lay)
     437             :       enddo
     438    12718080 :       do lay = laytrop+1, nlayers
     439   158699520 :          do ig = 1, ng2
     440   437944320 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + &
     441   583925760 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     442   145981440 :             taug(lay,ngs1+ig) = colh2o(lay) * &
     443   291962880 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     444   291962880 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     445   291962880 :                  fac01(lay) * absb(ind1(lay),ig) + &
     446   291962880 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
     447   729907200 :                  + taufor
     448   158146560 :             fracs(lay,ngs1+ig) = fracrefb(ig)
     449             :          enddo
     450             :       enddo
     451             : 
     452      552960 :       end subroutine taugb2
     453             : 
     454             : !----------------------------------------------------------------------------
     455      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers),  ig
     473     1105920 :       integer, dimension(nlayers) :: js, js1, jmn2o, jpl
     474     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
     475     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
     476     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2o, specparm_mn2o, specmult_mn2o, &
     477     1105920 :                        fmn2o, fmn2omf, chi_n2o, ratn2o, adjfac, adjcoln2o
     478     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
     479             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
     480     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
     481     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,9)/chi_mls(2,9)
     493             : 
     494             : !  P = 95.58 mb
     495      552960 :       refrat_planck_b = chi_mls(1,13)/chi_mls(2,13)
     496             : 
     497             : !  P = 706.270mb
     498      552960 :       refrat_m_a = chi_mls(1,3)/chi_mls(2,3)
     499             : 
     500             : !  P = 95.58 mb 
     501      552960 :       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    19906560 :       do lay = 1, laytrop
     510             : 
     511    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     512    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     513    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     514    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
     515    19353600 :          js(lay) = 1 + int(specmult(lay))
     516    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)        
     517             : 
     518    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     519    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     520    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     521    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
     522    19353600 :          js1(lay) = 1 + int(specmult1(lay))
     523    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     524             : 
     525    19353600 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
     526    19353600 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     527    19353600 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     528    19353600 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
     529    19353600 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     530    19353600 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     531    19353600 :          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    19353600 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     536    19906560 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     537             :       enddo
     538    19906560 :       do lay = 1, laytrop
     539    19906560 :          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    19353600 :             adjcoln2o(lay) = coln2o(lay)
     544             :          endif
     545             :       enddo
     546    19906560 :       do lay = 1, laytrop
     547             : 
     548    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     549    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     550    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     551    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     552    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
     553    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     554             : 
     555    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(3) + js(lay)
     556    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(3) + js1(lay)
     557             :       enddo
     558    19906560 :       do lay = 1, laytrop
     559             : 
     560    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
     561     4885892 :             p = fs(lay) - 1
     562     4885892 :             p4 = p**4
     563     4885892 :             fk0 = p4
     564     4885892 :             fk1 = 1 - p - 2.0_r8*p4
     565     4885892 :             fk2 = p + p4
     566     4885892 :             fac000(lay) = fk0*fac00(lay)
     567     4885892 :             fac100(lay) = fk1*fac00(lay)
     568     4885892 :             fac200(lay) = fk2*fac00(lay)
     569     4885892 :             fac010(lay) = fk0*fac10(lay)
     570     4885892 :             fac110(lay) = fk1*fac10(lay)
     571     4885892 :             fac210(lay) = fk2*fac10(lay)
     572    14467708 :          else if (specparm(lay) .gt. 0.875_r8) then
     573         735 :             p = -fs(lay) 
     574         735 :             p4 = p**4
     575         735 :             fk0 = p4
     576         735 :             fk1 = 1 - p - 2.0_r8*p4
     577         735 :             fk2 = p + p4
     578         735 :             fac000(lay) = fk0*fac00(lay)
     579         735 :             fac100(lay) = fk1*fac00(lay)
     580         735 :             fac200(lay) = fk2*fac00(lay)
     581         735 :             fac010(lay) = fk0*fac10(lay)
     582         735 :             fac110(lay) = fk1*fac10(lay)
     583         735 :             fac210(lay) = fk2*fac10(lay)
     584             :          else
     585    14466973 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     586    14466973 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     587    14466973 :             fac100(lay) = fs(lay) * fac00(lay)
     588    14466973 :             fac110(lay) = fs(lay) * fac10(lay)
     589             :          endif
     590    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
     591     2112458 :             p = fs1(lay) - 1
     592     2112458 :             p4 = p**4
     593     2112458 :             fk0 = p4
     594     2112458 :             fk1 = 1 - p - 2.0_r8*p4
     595     2112458 :             fk2 = p + p4
     596     2112458 :             fac001(lay) = fk0*fac01(lay)
     597     2112458 :             fac101(lay) = fk1*fac01(lay)
     598     2112458 :             fac201(lay) = fk2*fac01(lay)
     599     2112458 :             fac011(lay) = fk0*fac11(lay)
     600     2112458 :             fac111(lay) = fk1*fac11(lay)
     601     2112458 :             fac211(lay) = fk2*fac11(lay)
     602    17241142 :          else if (specparm1(lay) .gt. 0.875_r8) then
     603      171183 :             p = -fs1(lay) 
     604      171183 :             p4 = p**4
     605      171183 :             fk0 = p4
     606      171183 :             fk1 = 1 - p - 2.0_r8*p4
     607      171183 :             fk2 = p + p4
     608      171183 :             fac001(lay) = fk0*fac01(lay)
     609      171183 :             fac101(lay) = fk1*fac01(lay)
     610      171183 :             fac201(lay) = fk2*fac01(lay)
     611      171183 :             fac011(lay) = fk0*fac11(lay)
     612      171183 :             fac111(lay) = fk1*fac11(lay)
     613      171183 :             fac211(lay) = fk2*fac11(lay)
     614             :          else
     615    17069959 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     616    17069959 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     617    17069959 :             fac101(lay) = fs1(lay) * fac01(lay)
     618    17069959 :             fac111(lay) = fs1(lay) * fac11(lay)
     619             :          endif
     620             :       enddo
     621    19906560 :       do lay = 1, laytrop
     622             : 
     623   329564160 :          do ig = 1, ng3
     624   928972800 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     625  1238630400 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     626   619315200 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     627   619315200 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     628   928972800 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     629   928972800 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
     630   309657600 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     631   309657600 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     632   309657600 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     633             : 
     634   309657600 :             if (specparm(lay) .lt. 0.125_r8) then
     635             :                tau_major = speccomb(lay) * &
     636    78174272 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     637    78174272 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     638    78174272 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     639    78174272 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     640    78174272 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     641   469045632 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     642   231483328 :             else if (specparm(lay) .gt. 0.875_r8) then
     643             :                tau_major = speccomb(lay) * &
     644       11760 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     645       11760 :                     fac100(lay) * absa(ind0(lay),ig) + &
     646       11760 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     647       11760 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     648       11760 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     649       70560 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     650             :             else
     651             :                tau_major = speccomb(lay) * &
     652   231471568 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     653   231471568 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     654   231471568 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     655   925886272 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     656             :             endif
     657             : 
     658   309657600 :             if (specparm1(lay) .lt. 0.125_r8) then
     659             :                tau_major1 = speccomb1(lay) * &
     660    33799328 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     661    33799328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     662    33799328 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     663    33799328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     664    33799328 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     665   202795968 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     666   275858272 :             else if (specparm1(lay) .gt. 0.875_r8) then
     667             :                tau_major1 = speccomb1(lay) * &
     668     2738928 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     669     2738928 :                     fac101(lay) * absa(ind1(lay),ig) + &
     670     2738928 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     671     2738928 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     672     2738928 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     673    16433568 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     674             :             else
     675             :                tau_major1 = speccomb1(lay) * &
     676   273119344 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     677   273119344 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     678   273119344 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     679  1092477376 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     680             :             endif
     681             : 
     682   309657600 :             taug(lay,ngs2+ig) = tau_major + tau_major1 &
     683             :                  + tauself + taufor &
     684   309657600 :                  + adjcoln2o(lay)*absn2o
     685   619315200 :             fracs(lay,ngs2+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     686   638668800 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     687             :          enddo
     688             :       enddo
     689             : 
     690             : ! Upper atmosphere loop
     691    12718080 :       do lay = laytrop+1, nlayers
     692             : 
     693    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     694    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     695    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     696    12165120 :          specmult(lay) = 4._r8*(specparm(lay))
     697    12165120 :          js(lay) = 1 + int(specmult(lay))
     698    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
     699             : 
     700    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     701    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     702    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     703    12165120 :          specmult1(lay) = 4._r8*(specparm1(lay))
     704    12165120 :          js1(lay) = 1 + int(specmult1(lay))
     705    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     706             : 
     707    12165120 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     708    12165120 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     709    12165120 :          fac100(lay) = fs(lay) * fac00(lay)
     710    12165120 :          fac110(lay) = fs(lay) * fac10(lay)
     711    12165120 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     712    12165120 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     713    12165120 :          fac101(lay) = fs1(lay) * fac01(lay)
     714    12165120 :          fac111(lay) = fs1(lay) * fac11(lay)
     715             : 
     716    12165120 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_b*colco2(lay)
     717    12165120 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     718    12165120 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     719    12165120 :          specmult_mn2o(lay) = 4._r8*specparm_mn2o(lay)
     720    12165120 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     721    12165120 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     722    12165120 :          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    12165120 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     727    12718080 :          ratn2o(lay) = 1.e20*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     728             :       enddo
     729    12718080 :       do lay = laytrop+1, nlayers
     730    12718080 :          if (ratn2o(lay) .gt. 1.5_r8) then
     731     5638154 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
     732     5638154 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
     733             :          else
     734     6526966 :             adjcoln2o(lay) = coln2o(lay)
     735             :          endif
     736             :       enddo
     737    12718080 :       do lay = laytrop+1, nlayers
     738             : 
     739    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_b*colco2(lay)
     740    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     741    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     742    12165120 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
     743    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
     744    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     745             : 
     746    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(3) + js(lay)
     747    12718080 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(3) + js1(lay)
     748             :       enddo
     749    12718080 :       do lay = laytrop+1, nlayers
     750             : 
     751   207360000 :          do ig = 1, ng3
     752   583925760 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     753   778567680 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     754   583925760 :             n2om1 = kb_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     755   583925760 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay),ig)-kb_mn2o(jmn2o(lay),indminor(lay),ig))
     756   194641920 :             n2om2 = kb_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     757   194641920 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig)-kb_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     758   194641920 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     759   194641920 :             taug(lay,ngs2+ig) = speccomb(lay) * &
     760   194641920 :                 (fac000(lay) * absb(ind0(lay),ig) + &
     761   194641920 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
     762   194641920 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
     763   194641920 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
     764             :                 + speccomb1(lay) * &
     765   194641920 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
     766   194641920 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
     767   194641920 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
     768   194641920 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
     769             :                 + taufor &
     770  1751777280 :                 + adjcoln2o(lay)*absn2o
     771   389283840 :             fracs(lay,ngs2+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
     772   401448960 :                 (fracrefb(ig,jpl(lay)+1)-fracrefb(ig,jpl(lay)))
     773             :          enddo
     774             :       enddo
     775             : 
     776      552960 :       end subroutine taugb3
     777             : 
     778             : !----------------------------------------------------------------------------
     779      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
     796     1105920 :       integer, dimension(nlayers) :: js, js1, jpl
     797     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
     798     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
     799     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
     800             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
     801     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
     802     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,11)/chi_mls(2,11)
     810             : 
     811             : ! P = 95.58350 mb
     812      552960 :       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    19906560 :       do lay = 1, laytrop
     821             : 
     822    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     823    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     824    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     825    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
     826    19353600 :          js(lay) = 1 + int(specmult(lay))
     827    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
     828             : 
     829    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     830    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     831    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     832    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
     833    19353600 :          js1(lay) = 1 + int(specmult1(lay))
     834    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     835             : 
     836    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     837    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     838    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     839    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     840    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
     841    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     842             : 
     843    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(4) + js(lay)
     844    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(4) + js1(lay)
     845             :       enddo
     846    19906560 :       do lay = 1, laytrop
     847    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
     848     4885892 :             p = fs(lay) - 1
     849     4885892 :             p4 = p**4
     850     4885892 :             fk0 = p4
     851     4885892 :             fk1 = 1 - p - 2.0_r8*p4
     852     4885892 :             fk2 = p + p4
     853     4885892 :             fac000(lay) = fk0*fac00(lay)
     854     4885892 :             fac100(lay) = fk1*fac00(lay)
     855     4885892 :             fac200(lay) = fk2*fac00(lay)
     856     4885892 :             fac010(lay) = fk0*fac10(lay)
     857     4885892 :             fac110(lay) = fk1*fac10(lay)
     858     4885892 :             fac210(lay) = fk2*fac10(lay)
     859    14467708 :          else if (specparm(lay) .gt. 0.875_r8) then
     860         735 :             p = -fs(lay) 
     861         735 :             p4 = p**4
     862         735 :             fk0 = p4
     863         735 :             fk1 = 1 - p - 2.0_r8*p4
     864         735 :             fk2 = p + p4
     865         735 :             fac000(lay) = fk0*fac00(lay)
     866         735 :             fac100(lay) = fk1*fac00(lay)
     867         735 :             fac200(lay) = fk2*fac00(lay)
     868         735 :             fac010(lay) = fk0*fac10(lay)
     869         735 :             fac110(lay) = fk1*fac10(lay)
     870         735 :             fac210(lay) = fk2*fac10(lay)
     871             :          else
     872    14466973 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     873    14466973 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     874    14466973 :             fac100(lay) = fs(lay) * fac00(lay)
     875    14466973 :             fac110(lay) = fs(lay) * fac10(lay)
     876             :          endif
     877             : 
     878    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
     879     2112458 :             p = fs1(lay) - 1
     880     2112458 :             p4 = p**4
     881     2112458 :             fk0 = p4
     882     2112458 :             fk1 = 1 - p - 2.0_r8*p4
     883     2112458 :             fk2 = p + p4
     884     2112458 :             fac001(lay) = fk0*fac01(lay)
     885     2112458 :             fac101(lay) = fk1*fac01(lay)
     886     2112458 :             fac201(lay) = fk2*fac01(lay)
     887     2112458 :             fac011(lay) = fk0*fac11(lay)
     888     2112458 :             fac111(lay) = fk1*fac11(lay)
     889     2112458 :             fac211(lay) = fk2*fac11(lay)
     890    17241142 :          else if (specparm1(lay) .gt. 0.875_r8) then
     891      171183 :             p = -fs1(lay) 
     892      171183 :             p4 = p**4
     893      171183 :             fk0 = p4
     894      171183 :             fk1 = 1 - p - 2.0_r8*p4
     895      171183 :             fk2 = p + p4
     896      171183 :             fac001(lay) = fk0*fac01(lay)
     897      171183 :             fac101(lay) = fk1*fac01(lay)
     898      171183 :             fac201(lay) = fk2*fac01(lay)
     899      171183 :             fac011(lay) = fk0*fac11(lay)
     900      171183 :             fac111(lay) = fk1*fac11(lay)
     901      171183 :             fac211(lay) = fk2*fac11(lay)
     902             :          else
     903    17069959 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     904    17069959 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     905    17069959 :             fac101(lay) = fs1(lay) * fac01(lay)
     906    17069959 :             fac111(lay) = fs1(lay) * fac11(lay)
     907             :          endif
     908             :       enddo
     909    19906560 :       do lay = 1, laytrop
     910             : 
     911   290856960 :          do ig = 1, ng4
     912   812851200 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     913  1083801600 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     914   541900800 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     915   541900800 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     916             : 
     917   270950400 :             if (specparm(lay) .lt. 0.125_r8) then
     918             :                tau_major = speccomb(lay) * &
     919    68402488 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     920    68402488 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     921    68402488 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     922    68402488 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     923    68402488 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     924   410414928 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     925   202547912 :             else if (specparm(lay) .gt. 0.875_r8) then
     926             :                tau_major = speccomb(lay) * &
     927       10290 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     928       10290 :                     fac100(lay) * absa(ind0(lay),ig) + &
     929       10290 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     930       10290 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     931       10290 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     932       61740 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     933             :             else
     934             :                tau_major = speccomb(lay) * &
     935   202537622 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     936   202537622 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     937   202537622 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     938   810150488 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     939             :             endif
     940             : 
     941   270950400 :             if (specparm1(lay) .lt. 0.125_r8) then
     942             :                tau_major1 = speccomb1(lay) * &
     943    29574412 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     944    29574412 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     945    29574412 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     946    29574412 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     947    29574412 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     948   177446472 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     949   241375988 :             else if (specparm1(lay) .gt. 0.875_r8) then
     950             :                tau_major1 = speccomb1(lay) * &
     951     2396562 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     952     2396562 :                     fac101(lay) * absa(ind1(lay),ig) + &
     953     2396562 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     954     2396562 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     955     2396562 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     956    14379372 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     957             :             else
     958             :                tau_major1 = speccomb1(lay) * &
     959   238979426 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     960   238979426 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     961   238979426 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     962   955917704 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     963             :             endif
     964             : 
     965   270950400 :             taug(lay,ngs3+ig) = tau_major + tau_major1 &
     966   270950400 :                  + tauself + taufor
     967   541900800 :             fracs(lay,ngs3+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     968   561254400 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     969             :          enddo
     970             :       enddo
     971             : 
     972             : ! Upper atmosphere loop
     973    12718080 :       do lay = laytrop+1, nlayers
     974             : 
     975    12165120 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
     976    12165120 :          specparm(lay) = colo3(lay)/speccomb(lay)
     977    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     978    12165120 :          specmult(lay) = 4._r8*(specparm(lay))
     979    12165120 :          js(lay) = 1 + int(specmult(lay))
     980    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
     981             : 
     982    12165120 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
     983    12165120 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
     984    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     985    12165120 :          specmult1(lay) = 4._r8*(specparm1(lay))
     986    12165120 :          js1(lay) = 1 + int(specmult1(lay))
     987    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     988             : 
     989    12165120 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     990    12165120 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     991    12165120 :          fac100(lay) = fs(lay) * fac00(lay)
     992    12165120 :          fac110(lay) = fs(lay) * fac10(lay)
     993    12165120 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     994    12165120 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     995    12165120 :          fac101(lay) = fs1(lay) * fac01(lay)
     996    12165120 :          fac111(lay) = fs1(lay) * fac11(lay)
     997             : 
     998    12165120 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
     999    12165120 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1000    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1001    12165120 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1002    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1003    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1004             : 
    1005    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(4) + js(lay)
    1006    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(4) + js1(lay)
    1007   182476800 :          do ig = 1, ng4
    1008   170311680 :             taug(lay,ngs3+ig) =  speccomb(lay) * &
    1009   340623360 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1010   170311680 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1011   170311680 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1012   170311680 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1013             :                 + speccomb1(lay) * &
    1014   170311680 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
    1015   170311680 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1016   170311680 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1017  1703116800 :                 fac111(lay) * absb(ind1(lay)+6,ig))
    1018   340623360 :             fracs(lay,ngs3+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1019   352788480 :                 (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    12165120 :          taug(lay,ngs3+8)=taug(lay,ngs3+8)*0.92
    1026    12165120 :          taug(lay,ngs3+9)=taug(lay,ngs3+9)*0.88
    1027    12165120 :          taug(lay,ngs3+10)=taug(lay,ngs3+10)*1.07
    1028    12165120 :          taug(lay,ngs3+11)=taug(lay,ngs3+11)*1.1
    1029    12165120 :          taug(lay,ngs3+12)=taug(lay,ngs3+12)*0.99
    1030    12165120 :          taug(lay,ngs3+13)=taug(lay,ngs3+13)*0.88
    1031    12718080 :          taug(lay,ngs3+14)=taug(lay,ngs3+14)*0.943
    1032             : 
    1033             :       enddo
    1034             : 
    1035      552960 :       end subroutine taugb4
    1036             : 
    1037             : !----------------------------------------------------------------------------
    1038      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers),  ig
    1056     1105920 :       integer, dimension(nlayers) :: js, js1, jmo3, jpl
    1057     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1058     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1059     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mo3, specparm_mo3, specmult_mo3, fmo3
    1060     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1061             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1062     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1063     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,5)/chi_mls(2,5)
    1078             : 
    1079             : ! P = 0.2369 mb
    1080      552960 :       refrat_planck_b = chi_mls(3,43)/chi_mls(2,43)
    1081             : 
    1082             : ! P = 317.3480
    1083      552960 :       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    19906560 :       do lay = 1, laytrop
    1092             : 
    1093    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    1094    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1095    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1096    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    1097    19353600 :          js(lay) = 1 + int(specmult(lay))
    1098    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1099             : 
    1100    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    1101    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1102    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1103    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1104    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    1105    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1106             : 
    1107    19353600 :          speccomb_mo3(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
    1108    19353600 :          specparm_mo3(lay) = colh2o(lay)/speccomb_mo3(lay)
    1109    19353600 :          if (specparm_mo3(lay) .ge. oneminus) specparm_mo3(lay) = oneminus
    1110    19353600 :          specmult_mo3(lay) = 8._r8*specparm_mo3(lay)
    1111    19353600 :          jmo3(lay) = 1 + int(specmult_mo3(lay))
    1112    19353600 :          fmo3(lay) = mod(specmult_mo3(lay),1.0_r8)
    1113             : 
    1114    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    1115    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1116    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1117    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1118    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1119    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1120             : 
    1121    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(5) + js(lay)
    1122    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(5) + js1(lay)
    1123             :       enddo
    1124    19906560 :       do lay = 1, laytrop
    1125    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
    1126     4885892 :             p = fs(lay) - 1
    1127     4885892 :             p4 = p**4
    1128     4885892 :             fk0 = p4
    1129     4885892 :             fk1 = 1 - p - 2.0_r8*p4
    1130     4885892 :             fk2 = p + p4
    1131     4885892 :             fac000(lay) = fk0*fac00(lay)
    1132     4885892 :             fac100(lay) = fk1*fac00(lay)
    1133     4885892 :             fac200(lay) = fk2*fac00(lay)
    1134     4885892 :             fac010(lay) = fk0*fac10(lay)
    1135     4885892 :             fac110(lay) = fk1*fac10(lay)
    1136     4885892 :             fac210(lay) = fk2*fac10(lay)
    1137    14467708 :          else if (specparm(lay) .gt. 0.875_r8) then
    1138         735 :             p = -fs(lay) 
    1139         735 :             p4 = p**4
    1140         735 :             fk0 = p4
    1141         735 :             fk1 = 1 - p - 2.0_r8*p4
    1142         735 :             fk2 = p + p4
    1143         735 :             fac000(lay) = fk0*fac00(lay)
    1144         735 :             fac100(lay) = fk1*fac00(lay)
    1145         735 :             fac200(lay) = fk2*fac00(lay)
    1146         735 :             fac010(lay) = fk0*fac10(lay)
    1147         735 :             fac110(lay) = fk1*fac10(lay)
    1148         735 :             fac210(lay) = fk2*fac10(lay)
    1149             :          else
    1150    14466973 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1151    14466973 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1152    14466973 :             fac100(lay) = fs(lay) * fac00(lay)
    1153    14466973 :             fac110(lay) = fs(lay) * fac10(lay)
    1154             :          endif
    1155             : 
    1156    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
    1157     2112458 :             p = fs1(lay) - 1
    1158     2112458 :             p4 = p**4
    1159     2112458 :             fk0 = p4
    1160     2112458 :             fk1 = 1 - p - 2.0_r8*p4
    1161     2112458 :             fk2 = p + p4
    1162     2112458 :             fac001(lay) = fk0*fac01(lay)
    1163     2112458 :             fac101(lay) = fk1*fac01(lay)
    1164     2112458 :             fac201(lay) = fk2*fac01(lay)
    1165     2112458 :             fac011(lay) = fk0*fac11(lay)
    1166     2112458 :             fac111(lay) = fk1*fac11(lay)
    1167     2112458 :             fac211(lay) = fk2*fac11(lay)
    1168    17241142 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1169      171183 :             p = -fs1(lay) 
    1170      171183 :             p4 = p**4
    1171      171183 :             fk0 = p4
    1172      171183 :             fk1 = 1 - p - 2.0_r8*p4
    1173      171183 :             fk2 = p + p4
    1174      171183 :             fac001(lay) = fk0*fac01(lay)
    1175      171183 :             fac101(lay) = fk1*fac01(lay)
    1176      171183 :             fac201(lay) = fk2*fac01(lay)
    1177      171183 :             fac011(lay) = fk0*fac11(lay)
    1178      171183 :             fac111(lay) = fk1*fac11(lay)
    1179      171183 :             fac211(lay) = fk2*fac11(lay)
    1180             :          else
    1181    17069959 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1182    17069959 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1183    17069959 :             fac101(lay) = fs1(lay) * fac01(lay)
    1184    17069959 :             fac111(lay) = fs1(lay) * fac11(lay)
    1185             :          endif
    1186             :       enddo
    1187    19906560 :       do lay = 1, laytrop
    1188             : 
    1189   329564160 :          do ig = 1, ng5
    1190   928972800 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1191  1238630400 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1192   619315200 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1193   619315200 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1194   928972800 :             o3m1 = ka_mo3(jmo3(lay),indminor(lay),ig) + fmo3(lay) * &
    1195   928972800 :                  (ka_mo3(jmo3(lay)+1,indminor(lay),ig)-ka_mo3(jmo3(lay),indminor(lay),ig))
    1196   309657600 :             o3m2 = ka_mo3(jmo3(lay),indminor(lay)+1,ig) + fmo3(lay) * &
    1197   309657600 :                  (ka_mo3(jmo3(lay)+1,indminor(lay)+1,ig)-ka_mo3(jmo3(lay),indminor(lay)+1,ig))
    1198   309657600 :             abso3 = o3m1 + minorfrac(lay)*(o3m2-o3m1)
    1199             : 
    1200   309657600 :             if (specparm(lay) .lt. 0.125_r8) then
    1201             :                tau_major = speccomb(lay) * &
    1202    78174272 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1203    78174272 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1204    78174272 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1205    78174272 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1206    78174272 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1207   469045632 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1208   231483328 :             else if (specparm(lay) .gt. 0.875_r8) then
    1209             :                tau_major = speccomb(lay) * &
    1210       11760 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1211       11760 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1212       11760 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1213       11760 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1214       11760 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1215       70560 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1216             :             else
    1217             :                tau_major = speccomb(lay) * &
    1218   231471568 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1219   231471568 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1220   231471568 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1221   925886272 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1222             :             endif
    1223             : 
    1224   309657600 :             if (specparm1(lay) .lt. 0.125_r8) then
    1225             :                tau_major1 = speccomb1(lay) * &
    1226    33799328 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1227    33799328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1228    33799328 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1229    33799328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1230    33799328 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1231   202795968 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1232   275858272 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1233             :                tau_major1 = speccomb1(lay) * & 
    1234     2738928 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1235     2738928 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1236     2738928 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1237     2738928 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1238     2738928 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1239    16433568 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1240             :             else
    1241             :                tau_major1 = speccomb1(lay) * &
    1242   273119344 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1243   273119344 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1244   273119344 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1245  1092477376 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1246             :             endif
    1247             : 
    1248   309657600 :             taug(lay,ngs4+ig) = tau_major + tau_major1 &
    1249             :                  + tauself + taufor &
    1250   309657600 :                  + abso3*colo3(lay) &
    1251   619315200 :                  + wx(1,lay) * ccl4(ig)
    1252   619315200 :             fracs(lay,ngs4+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1253   638668800 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1254             :          enddo
    1255             :       enddo
    1256             : 
    1257             : ! Upper atmosphere loop
    1258    12718080 :       do lay = laytrop+1, nlayers
    1259             : 
    1260    12165120 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
    1261    12165120 :          specparm(lay) = colo3(lay)/speccomb(lay)
    1262    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1263    12165120 :          specmult(lay) = 4._r8*(specparm(lay))
    1264    12165120 :          js(lay) = 1 + int(specmult(lay))
    1265    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1266             : 
    1267    12165120 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
    1268    12165120 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
    1269    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1270    12165120 :          specmult1(lay) = 4._r8*(specparm1(lay))
    1271    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    1272    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1273             : 
    1274    12165120 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1275    12165120 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1276    12165120 :          fac100(lay) = fs(lay) * fac00(lay)
    1277    12165120 :          fac110(lay) = fs(lay) * fac10(lay)
    1278    12165120 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1279    12165120 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1280    12165120 :          fac101(lay) = fs1(lay) * fac01(lay)
    1281    12165120 :          fac111(lay) = fs1(lay) * fac11(lay)
    1282             : 
    1283    12165120 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
    1284    12165120 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1285    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1286    12165120 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1287    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1288    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1289             : 
    1290    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(5) + js(lay)
    1291    12718080 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(5) + js1(lay)
    1292             :       enddo
    1293    12718080 :       do lay = laytrop+1, nlayers
    1294   207360000 :          do ig = 1, ng5
    1295   389283840 :             taug(lay,ngs4+ig) = speccomb(lay) * &
    1296   389283840 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1297   194641920 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1298   194641920 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1299   194641920 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1300             :                 + speccomb1(lay) * &
    1301   194641920 :                 (fac001(lay) * absb(ind1(lay),ig) + &
    1302   194641920 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1303   194641920 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1304   194641920 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
    1305  2335703040 :                 + wx(1,lay) * ccl4(ig)
    1306   389283840 :             fracs(lay,ngs4+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1307   401448960 :                 (fracrefb(ig,jpl(lay)+1)-fracrefb(ig,jpl(lay)))
    1308             :          enddo
    1309             :       enddo
    1310             : 
    1311      552960 :       end subroutine taugb5
    1312             : 
    1313             : !----------------------------------------------------------------------------
    1314      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1332     1105920 :       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    19906560 :       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    19353600 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1351    19906560 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1352             :       enddo
    1353    19906560 :       do lay = 1, laytrop
    1354    19906560 :          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    19353600 :             adjcolco2(lay) = colco2(lay)
    1359             :          endif
    1360             :       enddo
    1361    19906560 :       do lay = 1, laytrop
    1362             : 
    1363    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(6) + 1
    1364    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(6) + 1
    1365             :       enddo
    1366    19906560 :       do lay = 1, laytrop
    1367             : 
    1368   174735360 :          do ig = 1, ng6
    1369   464486400 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1370   619315200 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1371   309657600 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1372   309657600 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1373   309657600 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1374   309657600 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1375   154828800 :             taug(lay,ngs5+ig) = colh2o(lay) * &
    1376   309657600 :                 (fac00(lay) * absa(ind0(lay),ig) + &
    1377   309657600 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1378   309657600 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1379   309657600 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    1380             :                  + tauself + taufor &
    1381             :                  + adjcolco2(lay) * absco2 &
    1382   154828800 :                  + wx(2,lay) * cfc11adj(ig) &
    1383   774144000 :                  + wx(3,lay) * cfc12(ig)
    1384   174182400 :             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    12718080 :       do lay = laytrop+1, nlayers
    1391             : 
    1392   110039040 :          do ig = 1, ng6
    1393   194641920 :             taug(lay,ngs5+ig) = 0.0_r8 &
    1394   194641920 :                  + wx(2,lay) * cfc11adj(ig) &
    1395   291962880 :                  + wx(3,lay) * cfc12(ig)
    1396   109486080 :             fracs(lay,ngs5+ig) = fracrefa(ig)
    1397             :          enddo
    1398             :       enddo
    1399             : 
    1400      552960 :       end subroutine taugb6
    1401             : 
    1402             : !----------------------------------------------------------------------------
    1403      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1421     1105920 :       integer, dimension(nlayers) :: js, js1, jmco2, jpl
    1422     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1423     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1424     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mco2, specparm_mco2, specmult_mco2, fmco2
    1425     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1426             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1427     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1428     1105920 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    1429             :       real(kind=r8) :: tauself, taufor, co2m1, co2m2, absco2
    1430     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,3)/chi_mls(3,3)
    1444             : 
    1445             : ! P = 706.2720 mb
    1446      552960 :       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    19906560 :       do lay = 1, laytrop
    1455             : 
    1456    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2oo3(lay)*colo3(lay)
    1457    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1458    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1459    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    1460    19353600 :          js(lay) = 1 + int(specmult(lay))
    1461    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1462             : 
    1463    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2oo3_1(lay)*colo3(lay)
    1464    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1465    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1466    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1467    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    1468    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1469             : 
    1470    19353600 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*colo3(lay)
    1471    19353600 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    1472    19353600 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    1473    19353600 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    1474             : 
    1475    19353600 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    1476    19353600 :          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    19353600 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1482    19353600 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1483    19353600 :          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    19353600 :             adjcolco2(lay) = colco2(lay)
    1488             :          endif
    1489             : 
    1490    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colo3(lay)
    1491    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1492    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1493    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1494    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1495    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1496             : 
    1497    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(7) + js(lay)
    1498    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(7) + js1(lay)
    1499             : 
    1500    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
    1501     4542062 :             p = fs(lay) - 1
    1502     4542062 :             p4 = p**4
    1503     4542062 :             fk0 = p4
    1504     4542062 :             fk1 = 1 - p - 2.0_r8*p4
    1505     4542062 :             fk2 = p + p4
    1506     4542062 :             fac000(lay) = fk0*fac00(lay)
    1507     4542062 :             fac100(lay) = fk1*fac00(lay)
    1508     4542062 :             fac200(lay) = fk2*fac00(lay)
    1509     4542062 :             fac010(lay) = fk0*fac10(lay)
    1510     4542062 :             fac110(lay) = fk1*fac10(lay)
    1511     4542062 :             fac210(lay) = fk2*fac10(lay)
    1512    14811538 :          else if (specparm(lay) .gt. 0.875_r8) then
    1513      585430 :             p = -fs(lay) 
    1514      585430 :             p4 = p**4
    1515      585430 :             fk0 = p4
    1516      585430 :             fk1 = 1 - p - 2.0_r8*p4
    1517      585430 :             fk2 = p + p4
    1518      585430 :             fac000(lay) = fk0*fac00(lay)
    1519      585430 :             fac100(lay) = fk1*fac00(lay)
    1520      585430 :             fac200(lay) = fk2*fac00(lay)
    1521      585430 :             fac010(lay) = fk0*fac10(lay)
    1522      585430 :             fac110(lay) = fk1*fac10(lay)
    1523      585430 :             fac210(lay) = fk2*fac10(lay)
    1524             :          else
    1525    14226108 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1526    14226108 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1527    14226108 :             fac100(lay) = fs(lay) * fac00(lay)
    1528    14226108 :             fac110(lay) = fs(lay) * fac10(lay)
    1529             :          endif
    1530    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
    1531     4542062 :             p = fs1(lay) - 1
    1532     4542062 :             p4 = p**4
    1533     4542062 :             fk0 = p4
    1534     4542062 :             fk1 = 1 - p - 2.0_r8*p4
    1535     4542062 :             fk2 = p + p4
    1536     4542062 :             fac001(lay) = fk0*fac01(lay)
    1537     4542062 :             fac101(lay) = fk1*fac01(lay)
    1538     4542062 :             fac201(lay) = fk2*fac01(lay)
    1539     4542062 :             fac011(lay) = fk0*fac11(lay)
    1540     4542062 :             fac111(lay) = fk1*fac11(lay)
    1541     4542062 :             fac211(lay) = fk2*fac11(lay)
    1542    14811538 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1543     1885294 :             p = -fs1(lay) 
    1544     1885294 :             p4 = p**4
    1545     1885294 :             fk0 = p4
    1546     1885294 :             fk1 = 1 - p - 2.0_r8*p4
    1547     1885294 :             fk2 = p + p4
    1548     1885294 :             fac001(lay) = fk0*fac01(lay)
    1549     1885294 :             fac101(lay) = fk1*fac01(lay)
    1550     1885294 :             fac201(lay) = fk2*fac01(lay)
    1551     1885294 :             fac011(lay) = fk0*fac11(lay)
    1552     1885294 :             fac111(lay) = fk1*fac11(lay)
    1553     1885294 :             fac211(lay) = fk2*fac11(lay)
    1554             :          else
    1555    12926244 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1556    12926244 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1557    12926244 :             fac101(lay) = fs1(lay) * fac01(lay)
    1558    12926244 :             fac111(lay) = fs1(lay) * fac11(lay)
    1559             :          endif
    1560             : 
    1561   252149760 :          do ig = 1, ng7
    1562   696729600 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1563   696729600 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1564   464486400 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1565   464486400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1566   696729600 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    1567   696729600 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    1568   232243200 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    1569   232243200 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    1570   232243200 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    1571             : 
    1572   232243200 :             if (specparm(lay) .lt. 0.125_r8) then
    1573             :                tau_major = speccomb(lay) * &
    1574    54504744 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1575    54504744 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1576    54504744 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1577    54504744 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1578    54504744 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1579   327028464 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1580   177738456 :             else if (specparm(lay) .gt. 0.875_r8) then
    1581             :                tau_major = speccomb(lay) * &
    1582     7025160 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1583     7025160 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1584     7025160 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1585     7025160 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1586     7025160 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1587    42150960 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1588             :             else
    1589             :                tau_major = speccomb(lay) * &
    1590   170713296 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1591   170713296 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1592   170713296 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1593   682853184 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1594             :             endif
    1595             : 
    1596   232243200 :             if (specparm1(lay) .lt. 0.125_r8) then
    1597             :                tau_major1 = speccomb1(lay) * &
    1598    22753068 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1599    22753068 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1600    22753068 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1601    22753068 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1602    22753068 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1603   136518408 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1604   209490132 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1605             :                tau_major1 = speccomb1(lay) * &
    1606    22623528 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1607    22623528 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1608    22623528 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1609    22623528 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1610    22623528 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1611   135741168 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1612             :             else
    1613             :                tau_major1 = speccomb1(lay) * &
    1614   186866604 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
    1615   186866604 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1616   186866604 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1617   747466416 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1618             :             endif
    1619             : 
    1620   232243200 :             taug(lay,ngs6+ig) = tau_major + tau_major1 &
    1621             :                  + tauself + taufor &
    1622   232243200 :                  + adjcolco2(lay)*absco2
    1623   464486400 :             fracs(lay,ngs6+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1624   483840000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1625             :          enddo
    1626             :       enddo
    1627             : 
    1628             : ! Upper atmosphere loop
    1629    12718080 :       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    12165120 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1635    12165120 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1636    12165120 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    1641             :          endif
    1642             : 
    1643    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(7) + 1
    1644    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(7) + 1
    1645             : 
    1646   158146560 :          do ig = 1, ng7
    1647   291962880 :             absco2 = kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1648   437944320 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig))
    1649   145981440 :             taug(lay,ngs6+ig) = colo3(lay) * &
    1650   291962880 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1651   291962880 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1652   291962880 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1653   291962880 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1654   729907200 :                  + adjcolco2(lay) * absco2
    1655   158146560 :             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    12165120 :          taug(lay,ngs6+6)=taug(lay,ngs6+6)*0.92_r8
    1662    12165120 :          taug(lay,ngs6+7)=taug(lay,ngs6+7)*0.88_r8
    1663    12165120 :          taug(lay,ngs6+8)=taug(lay,ngs6+8)*1.07_r8
    1664    12165120 :          taug(lay,ngs6+9)=taug(lay,ngs6+9)*1.1_r8
    1665    12165120 :          taug(lay,ngs6+10)=taug(lay,ngs6+10)*0.99_r8
    1666    12718080 :          taug(lay,ngs6+11)=taug(lay,ngs6+11)*0.855_r8
    1667             : 
    1668             :       enddo
    1669             : 
    1670      552960 :       end subroutine taugb7
    1671             : 
    1672             : !----------------------------------------------------------------------------
    1673      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1692             :       real(kind=r8) :: tauself, taufor, absco2, abso3, absn2o
    1693     1105920 :       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    19906560 :       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    19353600 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1716    19353600 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1717    19353600 :          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    19353600 :             adjcolco2(lay) = colco2(lay)
    1722             :          endif
    1723             : 
    1724    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(8) + 1
    1725    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(8) + 1
    1726             : 
    1727   174735360 :          do ig = 1, ng8
    1728   464486400 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1729   464486400 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1730   309657600 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1731   309657600 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1732   309657600 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1733   309657600 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1734             :             abso3 =  (ka_mo3(indminor(lay),ig) + minorfrac(lay) * &
    1735   154828800 :                  (ka_mo3(indminor(lay)+1,ig) - ka_mo3(indminor(lay),ig)))
    1736             :             absn2o =  (ka_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1737   154828800 :                  (ka_mn2o(indminor(lay)+1,ig) - ka_mn2o(indminor(lay),ig)))
    1738   154828800 :             taug(lay,ngs7+ig) = colh2o(lay) * &
    1739   309657600 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    1740   309657600 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1741   309657600 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1742   309657600 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    1743             :                  + tauself + taufor &
    1744             :                  + adjcolco2(lay)*absco2 &
    1745   154828800 :                  + colo3(lay) * abso3 &
    1746   154828800 :                  + coln2o(lay) * absn2o &
    1747   154828800 :                  + wx(3,lay) * cfc12(ig) &
    1748   774144000 :                  + wx(4,lay) * cfc22adj(ig)
    1749   174182400 :             fracs(lay,ngs7+ig) = fracrefa(ig)
    1750             :          enddo
    1751             :       enddo
    1752             : 
    1753             : ! Upper atmosphere loop
    1754    12718080 :       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    12165120 :          chi_co2(lay) = colco2(lay)/coldry(lay)
    1760    12165120 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1761    12165120 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    1766             :          endif
    1767             : 
    1768    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(8) + 1
    1769    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(8) + 1
    1770             : 
    1771   110039040 :          do ig = 1, ng8
    1772   194641920 :             absco2 =  (kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1773   291962880 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig)))
    1774             :             absn2o =  (kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1775    97320960 :                  (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig)))
    1776    97320960 :             taug(lay,ngs7+ig) = colo3(lay) * &
    1777   194641920 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1778   194641920 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1779   194641920 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1780   194641920 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1781             :                  + adjcolco2(lay)*absco2 &
    1782    97320960 :                  + coln2o(lay)*absn2o & 
    1783    97320960 :                  + wx(3,lay) * cfc12(ig) &
    1784   486604800 :                  + wx(4,lay) * cfc22adj(ig)
    1785   109486080 :             fracs(lay,ngs7+ig) = fracrefb(ig)
    1786             :          enddo
    1787             :       enddo
    1788             : 
    1789      552960 :       end subroutine taugb8
    1790             : 
    1791             : !----------------------------------------------------------------------------
    1792      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    1810     1105920 :       integer, dimension(nlayers) :: js, js1, jmn2o, jpl
    1811     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    1812     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    1813     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2o, specparm_mn2o, specmult_mn2o, fmn2o
    1814     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    1815             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    1816     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    1817     1105920 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    1818             :       real(kind=r8) :: tauself, taufor, n2om1, n2om2, absn2o
    1819     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,9)/chi_mls(6,9)
    1833             : 
    1834             : ! P = 706.272 mb 
    1835      552960 :       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    19906560 :       do lay = 1, laytrop
    1844             : 
    1845    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    1846    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1847    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1848    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    1849    19353600 :          js(lay) = 1 + int(specmult(lay))
    1850    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1851             : 
    1852    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    1853    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1854    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1855    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1856    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    1857    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1858             : 
    1859    19353600 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colch4(lay)
    1860    19353600 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
    1861    19353600 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
    1862    19353600 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
    1863    19353600 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
    1864    19353600 :          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    19353600 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    1870    19353600 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    1871    19353600 :          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    19353600 :             adjcoln2o(lay) = coln2o(lay)
    1876             :          endif
    1877             : 
    1878    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    1879    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1880    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1881    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1882    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1883    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1884             : 
    1885    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(9) + js(lay)
    1886    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(9) + js1(lay)
    1887             : 
    1888    19353600 :          if (specparm(lay) .lt. 0.125_r8) then
    1889     4847614 :             p = fs(lay) - 1
    1890     4847614 :             p4 = p**4
    1891     4847614 :             fk0 = p4
    1892     4847614 :             fk1 = 1 - p - 2.0_r8*p4
    1893     4847614 :             fk2 = p + p4
    1894     4847614 :             fac000(lay) = fk0*fac00(lay)
    1895     4847614 :             fac100(lay) = fk1*fac00(lay)
    1896     4847614 :             fac200(lay) = fk2*fac00(lay)
    1897     4847614 :             fac010(lay) = fk0*fac10(lay)
    1898     4847614 :             fac110(lay) = fk1*fac10(lay)
    1899     4847614 :             fac210(lay) = fk2*fac10(lay)
    1900    14505986 :          else if (specparm(lay) .gt. 0.875_r8) then
    1901         263 :             p = -fs(lay) 
    1902         263 :             p4 = p**4
    1903         263 :             fk0 = p4
    1904         263 :             fk1 = 1 - p - 2.0_r8*p4
    1905         263 :             fk2 = p + p4
    1906         263 :             fac000(lay) = fk0*fac00(lay)
    1907         263 :             fac100(lay) = fk1*fac00(lay)
    1908         263 :             fac200(lay) = fk2*fac00(lay)
    1909         263 :             fac010(lay) = fk0*fac10(lay)
    1910         263 :             fac110(lay) = fk1*fac10(lay)
    1911         263 :             fac210(lay) = fk2*fac10(lay)
    1912             :          else
    1913    14505723 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1914    14505723 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1915    14505723 :             fac100(lay) = fs(lay) * fac00(lay)
    1916    14505723 :             fac110(lay) = fs(lay) * fac10(lay)
    1917             :          endif
    1918             : 
    1919    19353600 :          if (specparm1(lay) .lt. 0.125_r8) then
    1920     2155639 :             p = fs1(lay) - 1
    1921     2155639 :             p4 = p**4
    1922     2155639 :             fk0 = p4
    1923     2155639 :             fk1 = 1 - p - 2.0_r8*p4
    1924     2155639 :             fk2 = p + p4
    1925     2155639 :             fac001(lay) = fk0*fac01(lay)
    1926     2155639 :             fac101(lay) = fk1*fac01(lay)
    1927     2155639 :             fac201(lay) = fk2*fac01(lay)
    1928     2155639 :             fac011(lay) = fk0*fac11(lay)
    1929     2155639 :             fac111(lay) = fk1*fac11(lay)
    1930     2155639 :             fac211(lay) = fk2*fac11(lay)
    1931    17197961 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1932      138311 :             p = -fs1(lay) 
    1933      138311 :             p4 = p**4
    1934      138311 :             fk0 = p4
    1935      138311 :             fk1 = 1 - p - 2.0_r8*p4
    1936      138311 :             fk2 = p + p4
    1937      138311 :             fac001(lay) = fk0*fac01(lay)
    1938      138311 :             fac101(lay) = fk1*fac01(lay)
    1939      138311 :             fac201(lay) = fk2*fac01(lay)
    1940      138311 :             fac011(lay) = fk0*fac11(lay)
    1941      138311 :             fac111(lay) = fk1*fac11(lay)
    1942      138311 :             fac211(lay) = fk2*fac11(lay)
    1943             :          else
    1944    17059650 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1945    17059650 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1946    17059650 :             fac101(lay) = fs1(lay) * fac01(lay)
    1947    17059650 :             fac111(lay) = fs1(lay) * fac11(lay)
    1948             :          endif
    1949             : 
    1950   252149760 :          do ig = 1, ng9
    1951   696729600 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1952   696729600 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1953   464486400 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1954   464486400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1955   696729600 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
    1956   696729600 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
    1957   232243200 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
    1958   232243200 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
    1959   232243200 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
    1960             : 
    1961   232243200 :             if (specparm(lay) .lt. 0.125_r8) then
    1962             :                tau_major = speccomb(lay) * &
    1963    58171368 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1964    58171368 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1965    58171368 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1966    58171368 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1967    58171368 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1968   349028208 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1969   174071832 :             else if (specparm(lay) .gt. 0.875_r8) then
    1970             :                tau_major = speccomb(lay) * &
    1971        3156 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1972        3156 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1973        3156 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1974        3156 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1975        3156 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1976       18936 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1977             :             else
    1978             :                tau_major = speccomb(lay) * &
    1979   174068676 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1980   174068676 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1981   174068676 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1982   696274704 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1983             :             endif
    1984             : 
    1985   232243200 :             if (specparm1(lay) .lt. 0.125_r8) then
    1986             :                tau_major1 = speccomb1(lay) * &
    1987    25867668 :                     (fac001(lay) * absa(ind1(lay),ig) + & 
    1988    25867668 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1989    25867668 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1990    25867668 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1991    25867668 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1992   155206008 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1993   206375532 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1994             :                tau_major1 = speccomb1(lay) * &
    1995     1659732 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1996     1659732 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1997     1659732 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1998     1659732 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1999     1659732 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2000     9958392 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2001             :             else
    2002             :                tau_major1 = speccomb1(lay) * &
    2003   204715800 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2004   204715800 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2005   204715800 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2006   818863200 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2007             :             endif
    2008             : 
    2009   232243200 :             taug(lay,ngs8+ig) = tau_major + tau_major1 &
    2010             :                  + tauself + taufor &
    2011   232243200 :                  + adjcoln2o(lay)*absn2o
    2012   464486400 :             fracs(lay,ngs8+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2013   483840000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2014             :          enddo
    2015             :       enddo
    2016             : 
    2017             : ! Upper atmosphere loop
    2018    12718080 :       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    12165120 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    2024    12165120 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    2025    12165120 :          if (ratn2o(lay) .gt. 1.5_r8) then
    2026     5638153 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
    2027     5638153 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
    2028             :          else
    2029     6526967 :             adjcoln2o(lay) = coln2o(lay)
    2030             :          endif
    2031             : 
    2032    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(9) + 1
    2033    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(9) + 1
    2034             : 
    2035   158699520 :          do ig = 1, ng9
    2036   291962880 :             absn2o = kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    2037   437944320 :                 (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig))
    2038   145981440 :             taug(lay,ngs8+ig) = colch4(lay) * &
    2039   291962880 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2040   291962880 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2041   291962880 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2042   291962880 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2043   729907200 :                  + adjcoln2o(lay)*absn2o
    2044   158146560 :             fracs(lay,ngs8+ig) = fracrefb(ig)
    2045             :          enddo
    2046             :       enddo
    2047             : 
    2048      552960 :       end subroutine taugb9
    2049             : 
    2050             : !----------------------------------------------------------------------------
    2051      552960 :       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     1105920 :       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    19906560 :       do lay = 1, laytrop
    2076    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(10) + 1
    2077    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(10) + 1
    2078             : 
    2079   136028160 :          do ig = 1, ng10
    2080   348364800 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2081   348364800 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2082   232243200 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2083   232243200 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2084   116121600 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2085   232243200 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2086   232243200 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2087   232243200 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2088   232243200 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    2089   580608000 :                  + tauself + taufor
    2090   135475200 :             fracs(lay,ngs9+ig) = fracrefa(ig)
    2091             :          enddo
    2092             :       enddo
    2093             : 
    2094             : ! Upper atmosphere loop
    2095    12718080 :       do lay = laytrop+1, nlayers
    2096    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(10) + 1
    2097    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(10) + 1
    2098             : 
    2099    85708800 :          do ig = 1, ng10
    2100   218972160 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2101   218972160 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2102    72990720 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2103   145981440 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2104   145981440 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2105   145981440 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2106   145981440 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2107   364953600 :                  + taufor
    2108    85155840 :             fracs(lay,ngs9+ig) = fracrefb(ig)
    2109             :          enddo
    2110             :       enddo
    2111             : 
    2112      552960 :       end subroutine taugb10
    2113             : 
    2114             : !----------------------------------------------------------------------------
    2115      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2132     1105920 :       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    19906560 :       do lay = 1, laytrop
    2145    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(11) + 1
    2146    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(11) + 1
    2147    19353600 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2148   174735360 :          do ig = 1, ng11
    2149   464486400 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2150   464486400 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2151   309657600 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2152   309657600 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    2153   309657600 :             tauo2 =  scaleo2(lay) * (ka_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2154   309657600 :                  (ka_mo2(indminor(lay)+1,ig) - ka_mo2(indminor(lay),ig)))
    2155   154828800 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2156   309657600 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2157   309657600 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2158   309657600 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2159   309657600 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2160             :                  + tauself + taufor &
    2161   774144000 :                  + tauo2
    2162   174182400 :             fracs(lay,ngs10+ig) = fracrefa(ig)
    2163             :          enddo
    2164             :       enddo
    2165             : 
    2166             : ! Upper atmosphere loop
    2167    12718080 :       do lay = laytrop+1, nlayers
    2168    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(11) + 1
    2169    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(11) + 1
    2170    12165120 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2171   110039040 :          do ig = 1, ng11
    2172   291962880 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2173   291962880 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2174   194641920 :             tauo2 =  scaleo2(lay) * (kb_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2175   194641920 :                  (kb_mo2(indminor(lay)+1,ig) - kb_mo2(indminor(lay),ig)))
    2176    97320960 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2177   194641920 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2178   194641920 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2179   194641920 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2180   194641920 :                  fac11(lay) * absb(ind1(lay)+1,ig))  &
    2181             :                  + taufor &
    2182   486604800 :                  + tauo2
    2183   109486080 :             fracs(lay,ngs10+ig) = fracrefb(ig)
    2184             :          enddo
    2185             :       enddo
    2186             : 
    2187      552960 :       end subroutine taugb11
    2188             : 
    2189             : !----------------------------------------------------------------------------
    2190      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2207     1105920 :       integer, dimension(nlayers) :: js, js1, jpl
    2208     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2209     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2210     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2211             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2212     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2213     1105920 :       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      552960 :       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    19906560 :       do lay = 1, laytrop
    2232             : 
    2233    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    2234    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2235    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2236    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    2237    19353600 :          js(lay) = 1 + int(specmult(lay))
    2238    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2239             : 
    2240    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    2241    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2242    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2243    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2244    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    2245    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2246             : 
    2247    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    2248    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2249    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2250    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2251    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2252    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2253             : 
    2254    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(12) + js(lay)
    2255    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(12) + js1(lay)
    2256             :       enddo
    2257    19906560 :       do lay = 1, laytrop
    2258    19906560 :          if (specparm(lay) .lt. 0.125_r8) then
    2259     4885892 :             p = fs(lay) - 1
    2260     4885892 :             p4 = p**4
    2261     4885892 :             fk0 = p4
    2262     4885892 :             fk1 = 1 - p - 2.0_r8*p4
    2263     4885892 :             fk2 = p + p4
    2264     4885892 :             fac000(lay) = fk0*fac00(lay)
    2265     4885892 :             fac100(lay) = fk1*fac00(lay)
    2266     4885892 :             fac200(lay) = fk2*fac00(lay)
    2267     4885892 :             fac010(lay) = fk0*fac10(lay)
    2268     4885892 :             fac110(lay) = fk1*fac10(lay)
    2269     4885892 :             fac210(lay) = fk2*fac10(lay)
    2270    14467708 :          else if (specparm(lay) .gt. 0.875_r8) then
    2271         735 :             p = -fs(lay) 
    2272         735 :             p4 = p**4
    2273         735 :             fk0 = p4
    2274         735 :             fk1 = 1 - p - 2.0_r8*p4
    2275         735 :             fk2 = p + p4
    2276         735 :             fac000(lay) = fk0*fac00(lay)
    2277         735 :             fac100(lay) = fk1*fac00(lay)
    2278         735 :             fac200(lay) = fk2*fac00(lay)
    2279         735 :             fac010(lay) = fk0*fac10(lay)
    2280         735 :             fac110(lay) = fk1*fac10(lay)
    2281         735 :             fac210(lay) = fk2*fac10(lay)
    2282             :          else
    2283    14466973 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2284    14466973 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2285    14466973 :             fac100(lay) = fs(lay) * fac00(lay)
    2286    14466973 :             fac110(lay) = fs(lay) * fac10(lay)
    2287             :          endif
    2288             : 
    2289             :       enddo
    2290    19906560 :       do lay = 1, laytrop
    2291    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
    2292     2112458 :             p = fs1(lay) - 1
    2293     2112458 :             p4 = p**4
    2294     2112458 :             fk0 = p4
    2295     2112458 :             fk1 = 1 - p - 2.0_r8*p4
    2296     2112458 :             fk2 = p + p4
    2297     2112458 :             fac001(lay) = fk0*fac01(lay)
    2298     2112458 :             fac101(lay) = fk1*fac01(lay)
    2299     2112458 :             fac201(lay) = fk2*fac01(lay)
    2300     2112458 :             fac011(lay) = fk0*fac11(lay)
    2301     2112458 :             fac111(lay) = fk1*fac11(lay)
    2302     2112458 :             fac211(lay) = fk2*fac11(lay)
    2303    17241142 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2304      171183 :             p = -fs1(lay) 
    2305      171183 :             p4 = p**4
    2306      171183 :             fk0 = p4
    2307      171183 :             fk1 = 1 - p - 2.0_r8*p4
    2308      171183 :             fk2 = p + p4
    2309      171183 :             fac001(lay) = fk0*fac01(lay)
    2310      171183 :             fac101(lay) = fk1*fac01(lay)
    2311      171183 :             fac201(lay) = fk2*fac01(lay)
    2312      171183 :             fac011(lay) = fk0*fac11(lay)
    2313      171183 :             fac111(lay) = fk1*fac11(lay)
    2314      171183 :             fac211(lay) = fk2*fac11(lay)
    2315             :          else
    2316    17069959 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2317    17069959 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2318    17069959 :             fac101(lay) = fs1(lay) * fac01(lay)
    2319    17069959 :             fac111(lay) = fs1(lay) * fac11(lay)
    2320             :          endif
    2321             :       enddo
    2322    19906560 :       do lay = 1, laytrop
    2323             : 
    2324   174735360 :          do ig = 1, ng12
    2325   464486400 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2326   619315200 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2327   309657600 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2328   309657600 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2329             : 
    2330   154828800 :             if (specparm(lay) .lt. 0.125_r8) then
    2331             :                tau_major = speccomb(lay) * &
    2332    39087136 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2333    39087136 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2334    39087136 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2335    39087136 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2336    39087136 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2337   234522816 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2338   115741664 :             else if (specparm(lay) .gt. 0.875_r8) then
    2339             :                tau_major = speccomb(lay) * &
    2340        5880 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2341        5880 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2342        5880 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2343        5880 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2344        5880 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2345       35280 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2346             :             else
    2347             :                tau_major = speccomb(lay) * &
    2348   115735784 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2349   115735784 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2350   115735784 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2351   462943136 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2352             :             endif
    2353             : 
    2354   154828800 :             if (specparm1(lay) .lt. 0.125_r8) then
    2355             :                tau_major1 = speccomb1(lay) * &
    2356    16899664 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2357    16899664 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2358    16899664 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2359    16899664 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2360    16899664 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2361   101397984 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2362   137929136 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2363             :                tau_major1 = speccomb1(lay) * &
    2364     1369464 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2365     1369464 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2366     1369464 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2367     1369464 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2368     1369464 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2369     8216784 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2370             :             else
    2371             :                tau_major1 = speccomb1(lay) * &
    2372   136559672 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2373   136559672 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2374   136559672 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2375   546238688 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2376             :             endif
    2377             : 
    2378   154828800 :             taug(lay,ngs11+ig) = tau_major + tau_major1 &
    2379   154828800 :                  + tauself + taufor
    2380   309657600 :             fracs(lay,ngs11+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2381   329011200 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2382             :          enddo
    2383             :       enddo
    2384             : 
    2385             : ! Upper atmosphere loop
    2386    12718080 :       do lay = laytrop+1, nlayers
    2387   110039040 :          do ig = 1, ng12
    2388    97320960 :             taug(lay,ngs11+ig) = 0.0_r8
    2389   109486080 :             fracs(lay,ngs11+ig) = 0.0_r8
    2390             :          enddo
    2391             :       enddo
    2392             : 
    2393      552960 :       end subroutine taugb12
    2394             : 
    2395             : !----------------------------------------------------------------------------
    2396      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2413     1105920 :       integer, dimension(nlayers) :: js, js1, jmco2, jmco, jpl
    2414     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2415     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2416     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mco2, specparm_mco2, specmult_mco2, fmco2
    2417     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mco, specparm_mco, specmult_mco, fmco
    2418     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2419             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2420     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2421     1105920 :       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     1105920 :       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      552960 :       refrat_planck_a = chi_mls(1,5)/chi_mls(4,5)
    2439             : 
    2440             : ! P = 1053. (Level 1)
    2441      552960 :       refrat_m_a = chi_mls(1,1)/chi_mls(4,1)
    2442             : 
    2443             : ! P = 706. (Level 3)
    2444      552960 :       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    19906560 :       do lay = 1, laytrop
    2453             : 
    2454    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2on2o(lay)*coln2o(lay)
    2455    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2456    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2457    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    2458    19353600 :          js(lay) = 1 + int(specmult(lay))
    2459    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2460             : 
    2461    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2on2o_1(lay)*coln2o(lay)
    2462    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2463    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2464    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2465    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    2466    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2467             : 
    2468    19353600 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*coln2o(lay)
    2469    19353600 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    2470    19353600 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    2471    19353600 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    2472    19353600 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    2473    19353600 :          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    19353600 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    2479    19906560 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/3.55e-4_r8
    2480             :       enddo
    2481    19906560 :       do lay = 1, laytrop
    2482    19906560 :          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    19353600 :             adjcolco2(lay) = colco2(lay)
    2487             :          endif
    2488             : 
    2489             :       enddo
    2490    19906560 :       do lay = 1, laytrop
    2491    19353600 :          speccomb_mco(lay) = colh2o(lay) + refrat_m_a3*coln2o(lay)
    2492    19353600 :          specparm_mco(lay) = colh2o(lay)/speccomb_mco(lay)
    2493    19353600 :          if (specparm_mco(lay) .ge. oneminus) specparm_mco(lay) = oneminus
    2494    19353600 :          specmult_mco(lay) = 8._r8*specparm_mco(lay)
    2495    19353600 :          jmco(lay) = 1 + int(specmult_mco(lay))
    2496    19353600 :          fmco(lay) = mod(specmult_mco(lay),1.0_r8)
    2497             : 
    2498    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*coln2o(lay)
    2499    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2500    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2501    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2502    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2503    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2504             : 
    2505    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(13) + js(lay)
    2506    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(13) + js1(lay)
    2507             :       enddo
    2508    19906560 :       do lay = 1, laytrop
    2509             : 
    2510    19906560 :          if (specparm(lay) .lt. 0.125_r8) then
    2511     4541335 :             p = fs(lay) - 1
    2512     4541335 :             p4 = p**4
    2513     4541335 :             fk0 = p4
    2514     4541335 :             fk1 = 1 - p - 2.0_r8*p4
    2515     4541335 :             fk2 = p + p4
    2516     4541335 :             fac000(lay) = fk0*fac00(lay)
    2517     4541335 :             fac100(lay) = fk1*fac00(lay)
    2518     4541335 :             fac200(lay) = fk2*fac00(lay)
    2519     4541335 :             fac010(lay) = fk0*fac10(lay)
    2520     4541335 :             fac110(lay) = fk1*fac10(lay)
    2521     4541335 :             fac210(lay) = fk2*fac10(lay)
    2522    14812265 :          else if (specparm(lay) .gt. 0.875_r8) then
    2523         560 :             p = -fs(lay) 
    2524         560 :             p4 = p**4
    2525         560 :             fk0 = p4
    2526         560 :             fk1 = 1 - p - 2.0_r8*p4
    2527         560 :             fk2 = p + p4
    2528         560 :             fac000(lay) = fk0*fac00(lay)
    2529         560 :             fac100(lay) = fk1*fac00(lay)
    2530         560 :             fac200(lay) = fk2*fac00(lay)
    2531         560 :             fac010(lay) = fk0*fac10(lay)
    2532         560 :             fac110(lay) = fk1*fac10(lay)
    2533         560 :             fac210(lay) = fk2*fac10(lay)
    2534             :          else
    2535    14811705 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2536    14811705 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2537    14811705 :             fac100(lay) = fs(lay) * fac00(lay)
    2538    14811705 :             fac110(lay) = fs(lay) * fac10(lay)
    2539             :          endif
    2540             : 
    2541             :       enddo
    2542    19906560 :       do lay = 1, laytrop
    2543    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
    2544     1939161 :             p = fs1(lay) - 1
    2545     1939161 :             p4 = p**4
    2546     1939161 :             fk0 = p4
    2547     1939161 :             fk1 = 1 - p - 2.0_r8*p4
    2548     1939161 :             fk2 = p + p4
    2549     1939161 :             fac001(lay) = fk0*fac01(lay)
    2550     1939161 :             fac101(lay) = fk1*fac01(lay)
    2551     1939161 :             fac201(lay) = fk2*fac01(lay)
    2552     1939161 :             fac011(lay) = fk0*fac11(lay)
    2553     1939161 :             fac111(lay) = fk1*fac11(lay)
    2554     1939161 :             fac211(lay) = fk2*fac11(lay)
    2555    17414439 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2556      176910 :             p = -fs1(lay) 
    2557      176910 :             p4 = p**4
    2558      176910 :             fk0 = p4
    2559      176910 :             fk1 = 1 - p - 2.0_r8*p4
    2560      176910 :             fk2 = p + p4
    2561      176910 :             fac001(lay) = fk0*fac01(lay)
    2562      176910 :             fac101(lay) = fk1*fac01(lay)
    2563      176910 :             fac201(lay) = fk2*fac01(lay)
    2564      176910 :             fac011(lay) = fk0*fac11(lay)
    2565      176910 :             fac111(lay) = fk1*fac11(lay)
    2566      176910 :             fac211(lay) = fk2*fac11(lay)
    2567             :          else
    2568    17237529 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2569    17237529 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2570    17237529 :             fac101(lay) = fs1(lay) * fac01(lay)
    2571    17237529 :             fac111(lay) = fs1(lay) * fac11(lay)
    2572             :          endif
    2573             :       enddo
    2574    19906560 :       do lay = 1, laytrop
    2575             : 
    2576    97320960 :          do ig = 1, ng13
    2577   232243200 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2578   309657600 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2579   154828800 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2580   154828800 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2581   232243200 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    2582   232243200 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    2583    77414400 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    2584    77414400 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    2585    77414400 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    2586    77414400 :             com1 = ka_mco(jmco(lay),indminor(lay),ig) + fmco(lay) * &
    2587   154828800 :                  (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    77414400 :                  (ka_mco(jmco(lay)+1,indminor(lay)+1,ig) - ka_mco(jmco(lay),indminor(lay)+1,ig))
    2590    77414400 :             absco = com1 + minorfrac(lay) * (com2 - com1)
    2591             : 
    2592    77414400 :             if (specparm(lay) .lt. 0.125_r8) then
    2593             :                tau_major = speccomb(lay) * &
    2594    18165340 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2595    18165340 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2596    18165340 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2597    18165340 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2598    18165340 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2599   108992040 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2600    59249060 :             else if (specparm(lay) .gt. 0.875_r8) then
    2601             :                tau_major = speccomb(lay) * &
    2602        2240 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2603        2240 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2604        2240 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2605        2240 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2606        2240 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2607       13440 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2608             :             else
    2609             :                tau_major = speccomb(lay) * &
    2610    59246820 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2611    59246820 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2612    59246820 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2613   236987280 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2614             :             endif
    2615             : 
    2616    77414400 :             if (specparm1(lay) .lt. 0.125_r8) then
    2617             :                tau_major1 = speccomb1(lay) * &
    2618     7756644 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2619     7756644 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2620     7756644 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2621     7756644 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2622     7756644 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2623    46539864 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2624    69657756 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2625             :                tau_major1 = speccomb1(lay) * &
    2626      707640 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2627      707640 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2628      707640 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2629      707640 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2630      707640 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2631     4245840 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2632             :             else
    2633             :                tau_major1 = speccomb1(lay) * &
    2634    68950116 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2635    68950116 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2636    68950116 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2637   275800464 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2638             :             endif
    2639             : 
    2640    77414400 :             taug(lay,ngs12+ig) = tau_major + tau_major1 &
    2641             :                  + tauself + taufor &
    2642             :                  + adjcolco2(lay)*absco2 &
    2643   154828800 :                  + colco(lay)*absco
    2644   154828800 :             fracs(lay,ngs12+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2645   174182400 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2646             :          enddo
    2647             :       enddo
    2648             : 
    2649             : ! Upper atmosphere loop
    2650    12718080 :       do lay = laytrop+1, nlayers
    2651    61378560 :          do ig = 1, ng13
    2652   145981440 :             abso3 = kb_mo3(indminor(lay),ig) + minorfrac(lay) * &
    2653   194641920 :                  (kb_mo3(indminor(lay)+1,ig) - kb_mo3(indminor(lay),ig))
    2654    48660480 :             taug(lay,ngs12+ig) = colo3(lay)*abso3
    2655    60825600 :             fracs(lay,ngs12+ig) =  fracrefb(ig)
    2656             :          enddo
    2657             :       enddo
    2658             : 
    2659      552960 :       end subroutine taugb13
    2660             : 
    2661             : !----------------------------------------------------------------------------
    2662      552960 :       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     1105920 :       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    19906560 :       do lay = 1, laytrop
    2687    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(14) + 1
    2688    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(14) + 1
    2689    58613760 :          do ig = 1, ng14
    2690   116121600 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2691   116121600 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2692    77414400 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2693    77414400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2694    38707200 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2695    77414400 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2696    77414400 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2697    77414400 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2698    77414400 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2699   193536000 :                  + tauself + taufor
    2700    58060800 :             fracs(lay,ngs13+ig) = fracrefa(ig)
    2701             :          enddo
    2702             :       enddo
    2703             : 
    2704             : ! Upper atmosphere loop
    2705    12718080 :       do lay = laytrop+1, nlayers
    2706    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(14) + 1
    2707    12165120 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(14) + 1
    2708    37048320 :          do ig = 1, ng14
    2709    24330240 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2710    72990720 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2711    48660480 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2712    48660480 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2713   145981440 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    2714    36495360 :             fracs(lay,ngs13+ig) = fracrefb(ig)
    2715             :          enddo
    2716             :       enddo
    2717             : 
    2718      552960 :       end subroutine taugb14
    2719             : 
    2720             : !----------------------------------------------------------------------------
    2721      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2739     1105920 :       integer, dimension(nlayers) :: js, js1, jmn2, jpl
    2740     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2741     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2742     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_mn2, specparm_mn2, specmult_mn2, fmn2
    2743     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2744             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2745     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2746     1105920 :       real(kind=r8), dimension(nlayers) :: fac001, fac101, fac201, fac011, fac111, fac211
    2747     1105920 :       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      552960 :       refrat_planck_a = chi_mls(4,1)/chi_mls(2,1)
    2759             : 
    2760             : ! P = 1053.
    2761      552960 :       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    19906560 :       do lay = 1, laytrop
    2770             : 
    2771    19353600 :          speccomb(lay) = coln2o(lay) + rat_n2oco2(lay)*colco2(lay)
    2772    19353600 :          specparm(lay) = coln2o(lay)/speccomb(lay)
    2773    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2774    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    2775    19353600 :          js(lay) = 1 + int(specmult(lay))
    2776    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2777             : 
    2778    19353600 :          speccomb1(lay) = coln2o(lay) + rat_n2oco2_1(lay)*colco2(lay)
    2779    19353600 :          specparm1(lay) = coln2o(lay)/speccomb1(lay)
    2780    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2781    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2782    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    2783    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2784             : 
    2785    19353600 :          speccomb_mn2(lay) = coln2o(lay) + refrat_m_a*colco2(lay)
    2786    19353600 :          specparm_mn2(lay) = coln2o(lay)/speccomb_mn2(lay)
    2787    19353600 :          if (specparm_mn2(lay) .ge. oneminus) specparm_mn2(lay) = oneminus
    2788    19353600 :          specmult_mn2(lay) = 8._r8*specparm_mn2(lay)
    2789    19353600 :          jmn2(lay) = 1 + int(specmult_mn2(lay))
    2790    19353600 :          fmn2(lay) = mod(specmult_mn2(lay),1.0_r8)
    2791             : 
    2792    19353600 :          speccomb_planck(lay) = coln2o(lay)+refrat_planck_a*colco2(lay)
    2793    19353600 :          specparm_planck(lay) = coln2o(lay)/speccomb_planck(lay)
    2794    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2795    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2796    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2797    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2798             : 
    2799    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(15) + js(lay)
    2800    19353600 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(15) + js1(lay)
    2801             :          
    2802    19906560 :          scalen2(lay) = colbrd(lay)*scaleminor(lay)
    2803             :       enddo
    2804    19906560 :       do lay = 1, laytrop
    2805    19906560 :          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    19353600 :          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    19353600 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2831    19353600 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2832    19353600 :             fac100(lay) = fs(lay) * fac00(lay)
    2833    19353600 :             fac110(lay) = fs(lay) * fac10(lay)
    2834             :          endif
    2835             :       enddo
    2836    19906560 :       do lay = 1, laytrop
    2837    19906560 :          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    19353600 :          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    19353600 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2863    19353600 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2864    19353600 :             fac101(lay) = fs1(lay) * fac01(lay)
    2865    19353600 :             fac111(lay) = fs1(lay) * fac11(lay)
    2866             :          endif
    2867             : 
    2868             :       enddo
    2869    19906560 :       do lay = 1, laytrop
    2870    58613760 :          do ig = 1, ng15
    2871   116121600 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2872   154828800 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2873    77414400 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2874    77414400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2875   116121600 :             n2m1 = ka_mn2(jmn2(lay),indminor(lay),ig) + fmn2(lay) * &
    2876   116121600 :                  (ka_mn2(jmn2(lay)+1,indminor(lay),ig) - ka_mn2(jmn2(lay),indminor(lay),ig))
    2877    38707200 :             n2m2 = ka_mn2(jmn2(lay),indminor(lay)+1,ig) + fmn2(lay) * &
    2878    38707200 :                  (ka_mn2(jmn2(lay)+1,indminor(lay)+1,ig) - ka_mn2(jmn2(lay),indminor(lay)+1,ig))
    2879    38707200 :             taun2 = scalen2(lay) * (n2m1 + minorfrac(lay) * (n2m2 - n2m1))
    2880             : 
    2881    38707200 :             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    38707200 :             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    38707200 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2900    38707200 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2901    38707200 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2902   154828800 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2903             :             endif 
    2904             : 
    2905    38707200 :             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    38707200 :             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    38707200 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2924    38707200 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2925    38707200 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2926   154828800 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2927             :             endif
    2928             : 
    2929    38707200 :             taug(lay,ngs14+ig) = tau_major + tau_major1 &
    2930             :                  + tauself + taufor &
    2931    38707200 :                  + taun2
    2932    77414400 :             fracs(lay,ngs14+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2933    96768000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2934             :          enddo
    2935             :       enddo
    2936             : 
    2937             : ! Upper atmosphere loop
    2938    12718080 :       do lay = laytrop+1, nlayers
    2939    37048320 :          do ig = 1, ng15
    2940    24330240 :             taug(lay,ngs14+ig) = 0.0_r8
    2941    36495360 :             fracs(lay,ngs14+ig) = 0.0_r8
    2942             :          enddo
    2943             :       enddo
    2944             : 
    2945      552960 :       end subroutine taugb15
    2946             : 
    2947             : !----------------------------------------------------------------------------
    2948      552960 :       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     1105920 :       integer :: lay, ind0(nlayers), ind1(nlayers), ig
    2965     1105920 :       integer, dimension(nlayers) :: js, js1, jpl
    2966     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb, specparm, specmult, fs
    2967     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb1, specparm1, specmult1, fs1
    2968     1105920 :       real(kind=r8), dimension(nlayers) :: speccomb_planck, specparm_planck, specmult_planck, fpl
    2969             :       real(kind=r8) :: p, p4, fk0, fk1, fk2
    2970     1105920 :       real(kind=r8), dimension(nlayers) :: fac000, fac100, fac200, fac010, fac110, fac210
    2971     1105920 :       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      552960 :       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    19906560 :       do lay = 1, laytrop
    2990             : 
    2991    19353600 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    2992    19353600 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2993    19353600 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2994    19353600 :          specmult(lay) = 8._r8*(specparm(lay))
    2995    19353600 :          js(lay) = 1 + int(specmult(lay))
    2996    19353600 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2997             : 
    2998    19353600 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    2999    19353600 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    3000    19353600 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    3001    19353600 :          specmult1(lay) = 8._r8*(specparm1(lay))
    3002    19353600 :          js1(lay) = 1 + int(specmult1(lay))
    3003    19353600 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    3004             : 
    3005    19353600 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    3006    19353600 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    3007    19353600 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    3008    19353600 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    3009    19353600 :          jpl(lay)= 1 + int(specmult_planck(lay))
    3010    19353600 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    3011             : 
    3012    19353600 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(16) + js(lay)
    3013    19906560 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(16) + js1(lay)
    3014             :       enddo
    3015    19906560 :       do lay = 1, laytrop
    3016    19906560 :          if (specparm(lay) .lt. 0.125_r8) then
    3017     4847614 :             p = fs(lay) - 1
    3018     4847614 :             p4 = p**4
    3019     4847614 :             fk0 = p4
    3020     4847614 :             fk1 = 1 - p - 2.0_r8*p4
    3021     4847614 :             fk2 = p + p4
    3022     4847614 :             fac000(lay) = fk0*fac00(lay)
    3023     4847614 :             fac100(lay) = fk1*fac00(lay)
    3024     4847614 :             fac200(lay) = fk2*fac00(lay)
    3025     4847614 :             fac010(lay) = fk0*fac10(lay)
    3026     4847614 :             fac110(lay) = fk1*fac10(lay)
    3027     4847614 :             fac210(lay) = fk2*fac10(lay)
    3028    14505986 :          else if (specparm(lay) .gt. 0.875_r8) then
    3029         263 :             p = -fs(lay) 
    3030         263 :             p4 = p**4
    3031         263 :             fk0 = p4
    3032         263 :             fk1 = 1 - p - 2.0_r8*p4
    3033         263 :             fk2 = p + p4
    3034         263 :             fac000(lay) = fk0*fac00(lay)
    3035         263 :             fac100(lay) = fk1*fac00(lay)
    3036         263 :             fac200(lay) = fk2*fac00(lay)
    3037         263 :             fac010(lay) = fk0*fac10(lay)
    3038         263 :             fac110(lay) = fk1*fac10(lay)
    3039         263 :             fac210(lay) = fk2*fac10(lay)
    3040             :          else
    3041    14505723 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    3042    14505723 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    3043    14505723 :             fac100(lay) = fs(lay) * fac00(lay)
    3044    14505723 :             fac110(lay) = fs(lay) * fac10(lay)
    3045             :          endif
    3046             :       enddo
    3047    19906560 :       do lay = 1, laytrop
    3048             : 
    3049    19906560 :          if (specparm1(lay) .lt. 0.125_r8) then
    3050     2155639 :             p = fs1(lay) - 1
    3051     2155639 :             p4 = p**4
    3052     2155639 :             fk0 = p4
    3053     2155639 :             fk1 = 1 - p - 2.0_r8*p4
    3054     2155639 :             fk2 = p + p4
    3055     2155639 :             fac001(lay) = fk0*fac01(lay)
    3056     2155639 :             fac101(lay) = fk1*fac01(lay)
    3057     2155639 :             fac201(lay) = fk2*fac01(lay)
    3058     2155639 :             fac011(lay) = fk0*fac11(lay)
    3059     2155639 :             fac111(lay) = fk1*fac11(lay)
    3060     2155639 :             fac211(lay) = fk2*fac11(lay)
    3061    17197961 :          else if (specparm1(lay) .gt. 0.875_r8) then
    3062      138311 :             p = -fs1(lay) 
    3063      138311 :             p4 = p**4
    3064      138311 :             fk0 = p4
    3065      138311 :             fk1 = 1 - p - 2.0_r8*p4
    3066      138311 :             fk2 = p + p4
    3067      138311 :             fac001(lay) = fk0*fac01(lay)
    3068      138311 :             fac101(lay) = fk1*fac01(lay)
    3069      138311 :             fac201(lay) = fk2*fac01(lay)
    3070      138311 :             fac011(lay) = fk0*fac11(lay)
    3071      138311 :             fac111(lay) = fk1*fac11(lay)
    3072      138311 :             fac211(lay) = fk2*fac11(lay)
    3073             :          else
    3074    17059650 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    3075    17059650 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    3076    17059650 :             fac101(lay) = fs1(lay) * fac01(lay)
    3077    17059650 :             fac111(lay) = fs1(lay) * fac11(lay)
    3078             :          endif
    3079             :       enddo
    3080    19906560 :       do lay = 1, laytrop
    3081             : 
    3082    58613760 :          do ig = 1, ng16
    3083   116121600 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    3084   154828800 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    3085    77414400 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    3086    77414400 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    3087             : 
    3088    38707200 :             if (specparm(lay) .lt. 0.125_r8) then
    3089             :                tau_major = speccomb(lay) * &
    3090     9695228 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3091     9695228 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3092     9695228 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    3093     9695228 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3094     9695228 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    3095    58171368 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    3096    29011972 :             else if (specparm(lay) .gt. 0.875_r8) then
    3097             :                tau_major = speccomb(lay) * &
    3098         526 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    3099         526 :                     fac100(lay) * absa(ind0(lay),ig) + &
    3100         526 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    3101         526 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    3102         526 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    3103        3156 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    3104             :             else
    3105             :                tau_major = speccomb(lay) * &
    3106    29011446 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3107    29011446 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3108    29011446 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3109   116045784 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    3110             :             endif
    3111             : 
    3112    38707200 :             if (specparm1(lay) .lt. 0.125_r8) then
    3113             :                tau_major1 = speccomb1(lay) * &
    3114     4311278 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3115     4311278 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3116     4311278 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    3117     4311278 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3118     4311278 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    3119    25867668 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    3120    34395922 :             else if (specparm1(lay) .gt. 0.875_r8) then
    3121             :                tau_major1 = speccomb1(lay) * &
    3122      276622 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    3123      276622 :                     fac101(lay) * absa(ind1(lay),ig) + &
    3124      276622 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    3125      276622 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    3126      276622 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    3127     1659732 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    3128             :             else
    3129             :                tau_major1 = speccomb1(lay) * &
    3130    34119300 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3131    34119300 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3132    34119300 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3133   136477200 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    3134             :             endif
    3135             : 
    3136    38707200 :             taug(lay,ngs15+ig) = tau_major + tau_major1 &
    3137    38707200 :                  + tauself + taufor
    3138    77414400 :             fracs(lay,ngs15+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    3139    96768000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    3140             :          enddo
    3141             :       enddo
    3142             : 
    3143             : ! Upper atmosphere loop
    3144    12718080 :       do lay = laytrop+1, nlayers
    3145    12165120 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(16) + 1
    3146    12718080 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(16) + 1
    3147             :       enddo
    3148    12718080 :       do lay = laytrop+1, nlayers
    3149    37048320 :          do ig = 1, ng16
    3150    48660480 :             taug(lay,ngs15+ig) = colch4(lay) * &
    3151    72990720 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    3152    48660480 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    3153    48660480 :                  fac01(lay) * absb(ind1(lay),ig) + &
    3154   170311680 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    3155    36495360 :             fracs(lay,ngs15+ig) = fracrefb(ig)
    3156             :          enddo
    3157             :       enddo
    3158             : 
    3159      552960 :       end subroutine taugb16
    3160             : 
    3161             :       end subroutine taumol
    3162             : 
    3163             :       end module rrtmg_lw_taumol
    3164             : 

Generated by: LCOV version 1.14