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:23:43 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    12718080 :       do lay = 1, laytrop
     315             : 
     316    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(1) + 1
     317    12165120 :          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    12165120 :          corradj(lay) =  1.
     323    12718080 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     324             :       enddo
     325    12718080 :       do lay = 1, laytrop
     326    12718080 :          if (pavel(lay) .lt. 250._r8) then
     327     3363387 :             corradj(lay) = 1._r8 - 0.15_r8 * (250._r8-pavel(lay)) / 154.4_r8
     328             :          endif
     329             :       enddo
     330    12718080 :       do lay = 1, laytrop
     331   134369280 :          do ig = 1, ng1
     332   364953600 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     333   486604800 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     334   243302400 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     335   243302400 :                  (forref(indfor(lay)+1,ig) -  forref(indfor(lay),ig))) 
     336   243302400 :             taun2 = scalen2(lay)*(ka_mn2(indminor(lay),ig) + & 
     337   243302400 :                  minorfrac(lay) * (ka_mn2(indminor(lay)+1,ig) - ka_mn2(indminor(lay),ig)))
     338   121651200 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     339   243302400 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     340   243302400 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     341   243302400 :                  fac01(lay) * absa(ind1(lay),ig) + &
     342   243302400 :                  fac11(lay) * absa(ind1(lay)+1,ig)) & 
     343   486604800 :                  + tauself + taufor + taun2)
     344   133816320 :              fracs(lay,ig) = fracrefa(ig)
     345             :          enddo
     346             :       enddo
     347             : 
     348             : ! Upper atmosphere loop
     349    18800640 :       do lay = laytrop+1, nlayers
     350             : 
     351    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(1) + 1
     352    18247680 :          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    18247680 :          corradj(lay) =  1._r8 - 0.15_r8 * (pavel(lay) / 95.6_r8)
     357             : 
     358    18800640 :          scalen2(lay) = colbrd(lay) * scaleminorn2(lay)
     359             :       enddo
     360    18800640 :       do lay = laytrop+1, nlayers
     361   201277440 :          do ig = 1, ng1
     362   547430400 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     363   729907200 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     364   364953600 :             taun2 = scalen2(lay)*(kb_mn2(indminor(lay),ig) + & 
     365   364953600 :                  minorfrac(lay) * (kb_mn2(indminor(lay)+1,ig) - kb_mn2(indminor(lay),ig)))
     366   182476800 :             taug(lay,ig) = corradj(lay) * (colh2o(lay) * &
     367   364953600 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     368   364953600 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     369   364953600 :                  fac01(lay) * absb(ind1(lay),ig) + &
     370   364953600 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &  
     371   729907200 :                  + taufor + taun2)
     372   200724480 :             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    12718080 :       do lay = 1, laytrop
     407             : 
     408    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(2) + 1
     409    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(2) + 1
     410             :          !inds = indself(lay)
     411             :          !indf = indfor(lay)
     412             :          !pp = pavel(lay)
     413    12718080 :          corradj(lay) = 1._r8 - .05_r8 * (pavel(lay) - 100._r8) / 900._r8
     414             :       enddo
     415    12718080 :       do lay = 1, laytrop
     416   158699520 :          do ig = 1, ng2
     417   437944320 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
     418   583925760 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     419   291962880 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     420   291962880 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     421   145981440 :             taug(lay,ngs1+ig) = corradj(lay) * (colh2o(lay) * &
     422   291962880 :                 (fac00(lay) * absa(ind0(lay),ig) + &
     423   291962880 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
     424   291962880 :                  fac01(lay) * absa(ind1(lay),ig) + &
     425   291962880 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
     426   729907200 :                  + tauself + taufor)
     427   158146560 :             fracs(lay,ngs1+ig) = fracrefa(ig)
     428             :          enddo
     429             :       enddo
     430             : 
     431             : ! Upper atmosphere loop
     432    18800640 :       do lay = laytrop+1, nlayers
     433             : 
     434    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(2) + 1
     435    18800640 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(2) + 1
     436             :          !indf = indfor(lay)
     437             :       enddo
     438    18800640 :       do lay = laytrop+1, nlayers
     439   237772800 :          do ig = 1, ng2
     440   656916480 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + &
     441   875888640 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     442   218972160 :             taug(lay,ngs1+ig) = colh2o(lay) * &
     443   437944320 :                 (fac00(lay) * absb(ind0(lay),ig) + &
     444   437944320 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
     445   437944320 :                  fac01(lay) * absb(ind1(lay),ig) + &
     446   437944320 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
     447  1094860800 :                  + taufor
     448   237219840 :             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    12718080 :       do lay = 1, laytrop
     510             : 
     511    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     512    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     513    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     514    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
     515    12165120 :          js(lay) = 1 + int(specmult(lay))
     516    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)        
     517             : 
     518    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     519    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     520    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     521    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
     522    12165120 :          js1(lay) = 1 + int(specmult1(lay))
     523    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     524             : 
     525    12165120 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
     526    12165120 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     527    12165120 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     528    12165120 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
     529    12165120 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     530    12165120 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     531    12165120 :          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    12165120 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     536    12718080 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     537             :       enddo
     538    12718080 :       do lay = 1, laytrop
     539    12718080 :          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    12165120 :             adjcoln2o(lay) = coln2o(lay)
     544             :          endif
     545             :       enddo
     546    12718080 :       do lay = 1, laytrop
     547             : 
     548    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     549    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     550    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     551    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     552    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
     553    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     554             : 
     555    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(3) + js(lay)
     556    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(3) + js1(lay)
     557             :       enddo
     558    12718080 :       do lay = 1, laytrop
     559             : 
     560    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
     561     2698634 :             p = fs(lay) - 1
     562     2698634 :             p4 = p**4
     563     2698634 :             fk0 = p4
     564     2698634 :             fk1 = 1 - p - 2.0_r8*p4
     565     2698634 :             fk2 = p + p4
     566     2698634 :             fac000(lay) = fk0*fac00(lay)
     567     2698634 :             fac100(lay) = fk1*fac00(lay)
     568     2698634 :             fac200(lay) = fk2*fac00(lay)
     569     2698634 :             fac010(lay) = fk0*fac10(lay)
     570     2698634 :             fac110(lay) = fk1*fac10(lay)
     571     2698634 :             fac210(lay) = fk2*fac10(lay)
     572     9466486 :          else if (specparm(lay) .gt. 0.875_r8) then
     573       13654 :             p = -fs(lay) 
     574       13654 :             p4 = p**4
     575       13654 :             fk0 = p4
     576       13654 :             fk1 = 1 - p - 2.0_r8*p4
     577       13654 :             fk2 = p + p4
     578       13654 :             fac000(lay) = fk0*fac00(lay)
     579       13654 :             fac100(lay) = fk1*fac00(lay)
     580       13654 :             fac200(lay) = fk2*fac00(lay)
     581       13654 :             fac010(lay) = fk0*fac10(lay)
     582       13654 :             fac110(lay) = fk1*fac10(lay)
     583       13654 :             fac210(lay) = fk2*fac10(lay)
     584             :          else
     585     9452832 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     586     9452832 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     587     9452832 :             fac100(lay) = fs(lay) * fac00(lay)
     588     9452832 :             fac110(lay) = fs(lay) * fac10(lay)
     589             :          endif
     590    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
     591     1146527 :             p = fs1(lay) - 1
     592     1146527 :             p4 = p**4
     593     1146527 :             fk0 = p4
     594     1146527 :             fk1 = 1 - p - 2.0_r8*p4
     595     1146527 :             fk2 = p + p4
     596     1146527 :             fac001(lay) = fk0*fac01(lay)
     597     1146527 :             fac101(lay) = fk1*fac01(lay)
     598     1146527 :             fac201(lay) = fk2*fac01(lay)
     599     1146527 :             fac011(lay) = fk0*fac11(lay)
     600     1146527 :             fac111(lay) = fk1*fac11(lay)
     601     1146527 :             fac211(lay) = fk2*fac11(lay)
     602    11018593 :          else if (specparm1(lay) .gt. 0.875_r8) then
     603      255010 :             p = -fs1(lay) 
     604      255010 :             p4 = p**4
     605      255010 :             fk0 = p4
     606      255010 :             fk1 = 1 - p - 2.0_r8*p4
     607      255010 :             fk2 = p + p4
     608      255010 :             fac001(lay) = fk0*fac01(lay)
     609      255010 :             fac101(lay) = fk1*fac01(lay)
     610      255010 :             fac201(lay) = fk2*fac01(lay)
     611      255010 :             fac011(lay) = fk0*fac11(lay)
     612      255010 :             fac111(lay) = fk1*fac11(lay)
     613      255010 :             fac211(lay) = fk2*fac11(lay)
     614             :          else
     615    10763583 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     616    10763583 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     617    10763583 :             fac101(lay) = fs1(lay) * fac01(lay)
     618    10763583 :             fac111(lay) = fs1(lay) * fac11(lay)
     619             :          endif
     620             :       enddo
     621    12718080 :       do lay = 1, laytrop
     622             : 
     623   207360000 :          do ig = 1, ng3
     624   583925760 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     625   778567680 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     626   389283840 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     627   389283840 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     628   583925760 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     629   583925760 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
     630   194641920 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     631   194641920 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     632   194641920 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     633             : 
     634   194641920 :             if (specparm(lay) .lt. 0.125_r8) then
     635             :                tau_major = speccomb(lay) * &
     636    43178144 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     637    43178144 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     638    43178144 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     639    43178144 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     640    43178144 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     641   259068864 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     642   151463776 :             else if (specparm(lay) .gt. 0.875_r8) then
     643             :                tau_major = speccomb(lay) * &
     644      218464 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     645      218464 :                     fac100(lay) * absa(ind0(lay),ig) + &
     646      218464 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     647      218464 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     648      218464 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     649     1310784 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     650             :             else
     651             :                tau_major = speccomb(lay) * &
     652   151245312 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     653   151245312 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     654   151245312 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     655   604981248 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     656             :             endif
     657             : 
     658   194641920 :             if (specparm1(lay) .lt. 0.125_r8) then
     659             :                tau_major1 = speccomb1(lay) * &
     660    18344432 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     661    18344432 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     662    18344432 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     663    18344432 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     664    18344432 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     665   110066592 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     666   176297488 :             else if (specparm1(lay) .gt. 0.875_r8) then
     667             :                tau_major1 = speccomb1(lay) * &
     668     4080160 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     669     4080160 :                     fac101(lay) * absa(ind1(lay),ig) + &
     670     4080160 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     671     4080160 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     672     4080160 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     673    24480960 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     674             :             else
     675             :                tau_major1 = speccomb1(lay) * &
     676   172217328 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     677   172217328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     678   172217328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     679   688869312 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     680             :             endif
     681             : 
     682   194641920 :             taug(lay,ngs2+ig) = tau_major + tau_major1 &
     683             :                  + tauself + taufor &
     684   194641920 :                  + adjcoln2o(lay)*absn2o
     685   389283840 :             fracs(lay,ngs2+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     686   401448960 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     687             :          enddo
     688             :       enddo
     689             : 
     690             : ! Upper atmosphere loop
     691    18800640 :       do lay = laytrop+1, nlayers
     692             : 
     693    18247680 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     694    18247680 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     695    18247680 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     696    18247680 :          specmult(lay) = 4._r8*(specparm(lay))
     697    18247680 :          js(lay) = 1 + int(specmult(lay))
     698    18247680 :          fs(lay) = mod(specmult(lay),1.0_r8)
     699             : 
     700    18247680 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     701    18247680 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     702    18247680 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     703    18247680 :          specmult1(lay) = 4._r8*(specparm1(lay))
     704    18247680 :          js1(lay) = 1 + int(specmult1(lay))
     705    18247680 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     706             : 
     707    18247680 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     708    18247680 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     709    18247680 :          fac100(lay) = fs(lay) * fac00(lay)
     710    18247680 :          fac110(lay) = fs(lay) * fac10(lay)
     711    18247680 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     712    18247680 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     713    18247680 :          fac101(lay) = fs1(lay) * fac01(lay)
     714    18247680 :          fac111(lay) = fs1(lay) * fac11(lay)
     715             : 
     716    18247680 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_b*colco2(lay)
     717    18247680 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
     718    18247680 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
     719    18247680 :          specmult_mn2o(lay) = 4._r8*specparm_mn2o(lay)
     720    18247680 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
     721    18247680 :          fmn2o(lay) = mod(specmult_mn2o(lay),1.0_r8)
     722    18247680 :          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    18247680 :          chi_n2o(lay) = coln2o(lay)/coldry(lay)
     727    18800640 :          ratn2o(lay) = 1.e20*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
     728             :       enddo
     729    18800640 :       do lay = laytrop+1, nlayers
     730    18800640 :          if (ratn2o(lay) .gt. 1.5_r8) then
     731     4897892 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
     732     4897892 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
     733             :          else
     734    13349788 :             adjcoln2o(lay) = coln2o(lay)
     735             :          endif
     736             :       enddo
     737    18800640 :       do lay = laytrop+1, nlayers
     738             : 
     739    18247680 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_b*colco2(lay)
     740    18247680 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     741    18247680 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     742    18247680 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
     743    18247680 :          jpl(lay)= 1 + int(specmult_planck(lay))
     744    18247680 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     745             : 
     746    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(3) + js(lay)
     747    18800640 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(3) + js1(lay)
     748             :       enddo
     749    18800640 :       do lay = laytrop+1, nlayers
     750             : 
     751   310763520 :          do ig = 1, ng3
     752   875888640 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + &
     753  1167851520 :                  forfrac(lay) * (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     754   875888640 :             n2om1 = kb_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
     755   875888640 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay),ig)-kb_mn2o(jmn2o(lay),indminor(lay),ig))
     756   291962880 :             n2om2 = kb_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
     757   291962880 :                  (kb_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig)-kb_mn2o(jmn2o(lay),indminor(lay)+1,ig))
     758   291962880 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
     759   291962880 :             taug(lay,ngs2+ig) = speccomb(lay) * &
     760   291962880 :                 (fac000(lay) * absb(ind0(lay),ig) + &
     761   291962880 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
     762   291962880 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
     763   291962880 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
     764             :                 + speccomb1(lay) * &
     765   291962880 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
     766   291962880 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
     767   291962880 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
     768   291962880 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
     769             :                 + taufor &
     770  2627665920 :                 + adjcoln2o(lay)*absn2o
     771   583925760 :             fracs(lay,ngs2+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
     772   602173440 :                 (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    12718080 :       do lay = 1, laytrop
     821             : 
     822    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
     823    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
     824    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     825    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
     826    12165120 :          js(lay) = 1 + int(specmult(lay))
     827    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
     828             : 
     829    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
     830    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
     831    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     832    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
     833    12165120 :          js1(lay) = 1 + int(specmult1(lay))
     834    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     835             : 
     836    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
     837    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
     838    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
     839    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
     840    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
     841    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
     842             : 
     843    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(4) + js(lay)
     844    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(4) + js1(lay)
     845             :       enddo
     846    12718080 :       do lay = 1, laytrop
     847    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
     848     2698634 :             p = fs(lay) - 1
     849     2698634 :             p4 = p**4
     850     2698634 :             fk0 = p4
     851     2698634 :             fk1 = 1 - p - 2.0_r8*p4
     852     2698634 :             fk2 = p + p4
     853     2698634 :             fac000(lay) = fk0*fac00(lay)
     854     2698634 :             fac100(lay) = fk1*fac00(lay)
     855     2698634 :             fac200(lay) = fk2*fac00(lay)
     856     2698634 :             fac010(lay) = fk0*fac10(lay)
     857     2698634 :             fac110(lay) = fk1*fac10(lay)
     858     2698634 :             fac210(lay) = fk2*fac10(lay)
     859     9466486 :          else if (specparm(lay) .gt. 0.875_r8) then
     860       13654 :             p = -fs(lay) 
     861       13654 :             p4 = p**4
     862       13654 :             fk0 = p4
     863       13654 :             fk1 = 1 - p - 2.0_r8*p4
     864       13654 :             fk2 = p + p4
     865       13654 :             fac000(lay) = fk0*fac00(lay)
     866       13654 :             fac100(lay) = fk1*fac00(lay)
     867       13654 :             fac200(lay) = fk2*fac00(lay)
     868       13654 :             fac010(lay) = fk0*fac10(lay)
     869       13654 :             fac110(lay) = fk1*fac10(lay)
     870       13654 :             fac210(lay) = fk2*fac10(lay)
     871             :          else
     872     9452832 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     873     9452832 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     874     9452832 :             fac100(lay) = fs(lay) * fac00(lay)
     875     9452832 :             fac110(lay) = fs(lay) * fac10(lay)
     876             :          endif
     877             : 
     878    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
     879     1146527 :             p = fs1(lay) - 1
     880     1146527 :             p4 = p**4
     881     1146527 :             fk0 = p4
     882     1146527 :             fk1 = 1 - p - 2.0_r8*p4
     883     1146527 :             fk2 = p + p4
     884     1146527 :             fac001(lay) = fk0*fac01(lay)
     885     1146527 :             fac101(lay) = fk1*fac01(lay)
     886     1146527 :             fac201(lay) = fk2*fac01(lay)
     887     1146527 :             fac011(lay) = fk0*fac11(lay)
     888     1146527 :             fac111(lay) = fk1*fac11(lay)
     889     1146527 :             fac211(lay) = fk2*fac11(lay)
     890    11018593 :          else if (specparm1(lay) .gt. 0.875_r8) then
     891      255010 :             p = -fs1(lay) 
     892      255010 :             p4 = p**4
     893      255010 :             fk0 = p4
     894      255010 :             fk1 = 1 - p - 2.0_r8*p4
     895      255010 :             fk2 = p + p4
     896      255010 :             fac001(lay) = fk0*fac01(lay)
     897      255010 :             fac101(lay) = fk1*fac01(lay)
     898      255010 :             fac201(lay) = fk2*fac01(lay)
     899      255010 :             fac011(lay) = fk0*fac11(lay)
     900      255010 :             fac111(lay) = fk1*fac11(lay)
     901      255010 :             fac211(lay) = fk2*fac11(lay)
     902             :          else
     903    10763583 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     904    10763583 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     905    10763583 :             fac101(lay) = fs1(lay) * fac01(lay)
     906    10763583 :             fac111(lay) = fs1(lay) * fac11(lay)
     907             :          endif
     908             :       enddo
     909    12718080 :       do lay = 1, laytrop
     910             : 
     911   183029760 :          do ig = 1, ng4
     912   510935040 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
     913   681246720 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
     914   340623360 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
     915   340623360 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
     916             : 
     917   170311680 :             if (specparm(lay) .lt. 0.125_r8) then
     918             :                tau_major = speccomb(lay) * &
     919    37780876 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     920    37780876 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     921    37780876 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
     922    37780876 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     923    37780876 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
     924   226685256 :                     fac210(lay) * absa(ind0(lay)+11,ig))
     925   132530804 :             else if (specparm(lay) .gt. 0.875_r8) then
     926             :                tau_major = speccomb(lay) * &
     927      191156 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
     928      191156 :                     fac100(lay) * absa(ind0(lay),ig) + &
     929      191156 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
     930      191156 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
     931      191156 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
     932     1146936 :                     fac010(lay) * absa(ind0(lay)+10,ig))
     933             :             else
     934             :                tau_major = speccomb(lay) * &
     935   132339648 :                     (fac000(lay) * absa(ind0(lay),ig) + &
     936   132339648 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
     937   132339648 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
     938   529358592 :                     fac110(lay) * absa(ind0(lay)+10,ig))
     939             :             endif
     940             : 
     941   170311680 :             if (specparm1(lay) .lt. 0.125_r8) then
     942             :                tau_major1 = speccomb1(lay) * &
     943    16051378 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
     944    16051378 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     945    16051378 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
     946    16051378 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     947    16051378 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
     948    96308268 :                     fac211(lay) * absa(ind1(lay)+11,ig))
     949   154260302 :             else if (specparm1(lay) .gt. 0.875_r8) then
     950             :                tau_major1 = speccomb1(lay) * &
     951     3570140 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
     952     3570140 :                     fac101(lay) * absa(ind1(lay),ig) + &
     953     3570140 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
     954     3570140 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
     955     3570140 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
     956    21420840 :                     fac011(lay) * absa(ind1(lay)+10,ig))
     957             :             else
     958             :                tau_major1 = speccomb1(lay) * &
     959   150690162 :                     (fac001(lay) * absa(ind1(lay),ig) + &
     960   150690162 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
     961   150690162 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
     962   602760648 :                     fac111(lay) * absa(ind1(lay)+10,ig))
     963             :             endif
     964             : 
     965   170311680 :             taug(lay,ngs3+ig) = tau_major + tau_major1 &
     966   170311680 :                  + tauself + taufor
     967   340623360 :             fracs(lay,ngs3+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
     968   352788480 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
     969             :          enddo
     970             :       enddo
     971             : 
     972             : ! Upper atmosphere loop
     973    18800640 :       do lay = laytrop+1, nlayers
     974             : 
     975    18247680 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
     976    18247680 :          specparm(lay) = colo3(lay)/speccomb(lay)
     977    18247680 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
     978    18247680 :          specmult(lay) = 4._r8*(specparm(lay))
     979    18247680 :          js(lay) = 1 + int(specmult(lay))
     980    18247680 :          fs(lay) = mod(specmult(lay),1.0_r8)
     981             : 
     982    18247680 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
     983    18247680 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
     984    18247680 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
     985    18247680 :          specmult1(lay) = 4._r8*(specparm1(lay))
     986    18247680 :          js1(lay) = 1 + int(specmult1(lay))
     987    18247680 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
     988             : 
     989    18247680 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
     990    18247680 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
     991    18247680 :          fac100(lay) = fs(lay) * fac00(lay)
     992    18247680 :          fac110(lay) = fs(lay) * fac10(lay)
     993    18247680 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
     994    18247680 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
     995    18247680 :          fac101(lay) = fs1(lay) * fac01(lay)
     996    18247680 :          fac111(lay) = fs1(lay) * fac11(lay)
     997             : 
     998    18247680 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
     999    18247680 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1000    18247680 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1001    18247680 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1002    18247680 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1003    18247680 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1004             : 
    1005    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(4) + js(lay)
    1006    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(4) + js1(lay)
    1007   273715200 :          do ig = 1, ng4
    1008   255467520 :             taug(lay,ngs3+ig) =  speccomb(lay) * &
    1009   510935040 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1010   255467520 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1011   255467520 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1012   255467520 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1013             :                 + speccomb1(lay) * &
    1014   255467520 :                 (fac001(lay) * absb(ind1(lay),ig) +  &
    1015   255467520 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1016   255467520 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1017  2554675200 :                 fac111(lay) * absb(ind1(lay)+6,ig))
    1018   510935040 :             fracs(lay,ngs3+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1019   529182720 :                 (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    18247680 :          taug(lay,ngs3+8)=taug(lay,ngs3+8)*0.92
    1026    18247680 :          taug(lay,ngs3+9)=taug(lay,ngs3+9)*0.88
    1027    18247680 :          taug(lay,ngs3+10)=taug(lay,ngs3+10)*1.07
    1028    18247680 :          taug(lay,ngs3+11)=taug(lay,ngs3+11)*1.1
    1029    18247680 :          taug(lay,ngs3+12)=taug(lay,ngs3+12)*0.99
    1030    18247680 :          taug(lay,ngs3+13)=taug(lay,ngs3+13)*0.88
    1031    18800640 :          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    12718080 :       do lay = 1, laytrop
    1092             : 
    1093    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    1094    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1095    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1096    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    1097    12165120 :          js(lay) = 1 + int(specmult(lay))
    1098    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1099             : 
    1100    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    1101    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1102    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1103    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1104    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    1105    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1106             : 
    1107    12165120 :          speccomb_mo3(lay) = colh2o(lay) + refrat_m_a*colco2(lay)
    1108    12165120 :          specparm_mo3(lay) = colh2o(lay)/speccomb_mo3(lay)
    1109    12165120 :          if (specparm_mo3(lay) .ge. oneminus) specparm_mo3(lay) = oneminus
    1110    12165120 :          specmult_mo3(lay) = 8._r8*specparm_mo3(lay)
    1111    12165120 :          jmo3(lay) = 1 + int(specmult_mo3(lay))
    1112    12165120 :          fmo3(lay) = mod(specmult_mo3(lay),1.0_r8)
    1113             : 
    1114    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    1115    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1116    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1117    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1118    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1119    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1120             : 
    1121    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(5) + js(lay)
    1122    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(5) + js1(lay)
    1123             :       enddo
    1124    12718080 :       do lay = 1, laytrop
    1125    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
    1126     2698634 :             p = fs(lay) - 1
    1127     2698634 :             p4 = p**4
    1128     2698634 :             fk0 = p4
    1129     2698634 :             fk1 = 1 - p - 2.0_r8*p4
    1130     2698634 :             fk2 = p + p4
    1131     2698634 :             fac000(lay) = fk0*fac00(lay)
    1132     2698634 :             fac100(lay) = fk1*fac00(lay)
    1133     2698634 :             fac200(lay) = fk2*fac00(lay)
    1134     2698634 :             fac010(lay) = fk0*fac10(lay)
    1135     2698634 :             fac110(lay) = fk1*fac10(lay)
    1136     2698634 :             fac210(lay) = fk2*fac10(lay)
    1137     9466486 :          else if (specparm(lay) .gt. 0.875_r8) then
    1138       13654 :             p = -fs(lay) 
    1139       13654 :             p4 = p**4
    1140       13654 :             fk0 = p4
    1141       13654 :             fk1 = 1 - p - 2.0_r8*p4
    1142       13654 :             fk2 = p + p4
    1143       13654 :             fac000(lay) = fk0*fac00(lay)
    1144       13654 :             fac100(lay) = fk1*fac00(lay)
    1145       13654 :             fac200(lay) = fk2*fac00(lay)
    1146       13654 :             fac010(lay) = fk0*fac10(lay)
    1147       13654 :             fac110(lay) = fk1*fac10(lay)
    1148       13654 :             fac210(lay) = fk2*fac10(lay)
    1149             :          else
    1150     9452832 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1151     9452832 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1152     9452832 :             fac100(lay) = fs(lay) * fac00(lay)
    1153     9452832 :             fac110(lay) = fs(lay) * fac10(lay)
    1154             :          endif
    1155             : 
    1156    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
    1157     1146527 :             p = fs1(lay) - 1
    1158     1146527 :             p4 = p**4
    1159     1146527 :             fk0 = p4
    1160     1146527 :             fk1 = 1 - p - 2.0_r8*p4
    1161     1146527 :             fk2 = p + p4
    1162     1146527 :             fac001(lay) = fk0*fac01(lay)
    1163     1146527 :             fac101(lay) = fk1*fac01(lay)
    1164     1146527 :             fac201(lay) = fk2*fac01(lay)
    1165     1146527 :             fac011(lay) = fk0*fac11(lay)
    1166     1146527 :             fac111(lay) = fk1*fac11(lay)
    1167     1146527 :             fac211(lay) = fk2*fac11(lay)
    1168    11018593 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1169      255010 :             p = -fs1(lay) 
    1170      255010 :             p4 = p**4
    1171      255010 :             fk0 = p4
    1172      255010 :             fk1 = 1 - p - 2.0_r8*p4
    1173      255010 :             fk2 = p + p4
    1174      255010 :             fac001(lay) = fk0*fac01(lay)
    1175      255010 :             fac101(lay) = fk1*fac01(lay)
    1176      255010 :             fac201(lay) = fk2*fac01(lay)
    1177      255010 :             fac011(lay) = fk0*fac11(lay)
    1178      255010 :             fac111(lay) = fk1*fac11(lay)
    1179      255010 :             fac211(lay) = fk2*fac11(lay)
    1180             :          else
    1181    10763583 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1182    10763583 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1183    10763583 :             fac101(lay) = fs1(lay) * fac01(lay)
    1184    10763583 :             fac111(lay) = fs1(lay) * fac11(lay)
    1185             :          endif
    1186             :       enddo
    1187    12718080 :       do lay = 1, laytrop
    1188             : 
    1189   207360000 :          do ig = 1, ng5
    1190   583925760 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1191   778567680 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1192   389283840 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1193   389283840 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1194   583925760 :             o3m1 = ka_mo3(jmo3(lay),indminor(lay),ig) + fmo3(lay) * &
    1195   583925760 :                  (ka_mo3(jmo3(lay)+1,indminor(lay),ig)-ka_mo3(jmo3(lay),indminor(lay),ig))
    1196   194641920 :             o3m2 = ka_mo3(jmo3(lay),indminor(lay)+1,ig) + fmo3(lay) * &
    1197   194641920 :                  (ka_mo3(jmo3(lay)+1,indminor(lay)+1,ig)-ka_mo3(jmo3(lay),indminor(lay)+1,ig))
    1198   194641920 :             abso3 = o3m1 + minorfrac(lay)*(o3m2-o3m1)
    1199             : 
    1200   194641920 :             if (specparm(lay) .lt. 0.125_r8) then
    1201             :                tau_major = speccomb(lay) * &
    1202    43178144 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1203    43178144 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1204    43178144 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1205    43178144 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1206    43178144 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1207   259068864 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1208   151463776 :             else if (specparm(lay) .gt. 0.875_r8) then
    1209             :                tau_major = speccomb(lay) * &
    1210      218464 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1211      218464 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1212      218464 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1213      218464 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1214      218464 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1215     1310784 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1216             :             else
    1217             :                tau_major = speccomb(lay) * &
    1218   151245312 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1219   151245312 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1220   151245312 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1221   604981248 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1222             :             endif
    1223             : 
    1224   194641920 :             if (specparm1(lay) .lt. 0.125_r8) then
    1225             :                tau_major1 = speccomb1(lay) * &
    1226    18344432 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1227    18344432 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1228    18344432 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1229    18344432 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1230    18344432 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1231   110066592 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1232   176297488 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1233             :                tau_major1 = speccomb1(lay) * & 
    1234     4080160 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1235     4080160 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1236     4080160 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1237     4080160 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1238     4080160 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1239    24480960 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1240             :             else
    1241             :                tau_major1 = speccomb1(lay) * &
    1242   172217328 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1243   172217328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1244   172217328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1245   688869312 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1246             :             endif
    1247             : 
    1248   194641920 :             taug(lay,ngs4+ig) = tau_major + tau_major1 &
    1249             :                  + tauself + taufor &
    1250   194641920 :                  + abso3*colo3(lay) &
    1251   389283840 :                  + wx(1,lay) * ccl4(ig)
    1252   389283840 :             fracs(lay,ngs4+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1253   401448960 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1254             :          enddo
    1255             :       enddo
    1256             : 
    1257             : ! Upper atmosphere loop
    1258    18800640 :       do lay = laytrop+1, nlayers
    1259             : 
    1260    18247680 :          speccomb(lay) = colo3(lay) + rat_o3co2(lay)*colco2(lay)
    1261    18247680 :          specparm(lay) = colo3(lay)/speccomb(lay)
    1262    18247680 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1263    18247680 :          specmult(lay) = 4._r8*(specparm(lay))
    1264    18247680 :          js(lay) = 1 + int(specmult(lay))
    1265    18247680 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1266             : 
    1267    18247680 :          speccomb1(lay) = colo3(lay) + rat_o3co2_1(lay)*colco2(lay)
    1268    18247680 :          specparm1(lay) = colo3(lay)/speccomb1(lay)
    1269    18247680 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1270    18247680 :          specmult1(lay) = 4._r8*(specparm1(lay))
    1271    18247680 :          js1(lay) = 1 + int(specmult1(lay))
    1272    18247680 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1273             : 
    1274    18247680 :          fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1275    18247680 :          fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1276    18247680 :          fac100(lay) = fs(lay) * fac00(lay)
    1277    18247680 :          fac110(lay) = fs(lay) * fac10(lay)
    1278    18247680 :          fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1279    18247680 :          fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1280    18247680 :          fac101(lay) = fs1(lay) * fac01(lay)
    1281    18247680 :          fac111(lay) = fs1(lay) * fac11(lay)
    1282             : 
    1283    18247680 :          speccomb_planck(lay) = colo3(lay)+refrat_planck_b*colco2(lay)
    1284    18247680 :          specparm_planck(lay) = colo3(lay)/speccomb_planck(lay)
    1285    18247680 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1286    18247680 :          specmult_planck(lay) = 4._r8*specparm_planck(lay)
    1287    18247680 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1288    18247680 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1289             : 
    1290    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(5) + js(lay)
    1291    18800640 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(5) + js1(lay)
    1292             :       enddo
    1293    18800640 :       do lay = laytrop+1, nlayers
    1294   310763520 :          do ig = 1, ng5
    1295   583925760 :             taug(lay,ngs4+ig) = speccomb(lay) * &
    1296   583925760 :                 (fac000(lay) * absb(ind0(lay),ig) + &
    1297   291962880 :                 fac100(lay) * absb(ind0(lay)+1,ig) + &
    1298   291962880 :                 fac010(lay) * absb(ind0(lay)+5,ig) + &
    1299   291962880 :                 fac110(lay) * absb(ind0(lay)+6,ig)) &
    1300             :                 + speccomb1(lay) * &
    1301   291962880 :                 (fac001(lay) * absb(ind1(lay),ig) + &
    1302   291962880 :                 fac101(lay) * absb(ind1(lay)+1,ig) + &
    1303   291962880 :                 fac011(lay) * absb(ind1(lay)+5,ig) + &
    1304   291962880 :                 fac111(lay) * absb(ind1(lay)+6,ig))  &
    1305  3503554560 :                 + wx(1,lay) * ccl4(ig)
    1306   583925760 :             fracs(lay,ngs4+ig) = fracrefb(ig,jpl(lay)) + fpl(lay) * &
    1307   602173440 :                 (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    12718080 :       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    12165120 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1351    12718080 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1352             :       enddo
    1353    12718080 :       do lay = 1, laytrop
    1354    12718080 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    1359             :          endif
    1360             :       enddo
    1361    12718080 :       do lay = 1, laytrop
    1362             : 
    1363    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(6) + 1
    1364    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(6) + 1
    1365             :       enddo
    1366    12718080 :       do lay = 1, laytrop
    1367             : 
    1368   110039040 :          do ig = 1, ng6
    1369   291962880 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1370   389283840 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1371   194641920 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1372   194641920 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1373   194641920 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1374   194641920 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1375    97320960 :             taug(lay,ngs5+ig) = colh2o(lay) * &
    1376   194641920 :                 (fac00(lay) * absa(ind0(lay),ig) + &
    1377   194641920 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1378   194641920 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1379   194641920 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    1380             :                  + tauself + taufor &
    1381             :                  + adjcolco2(lay) * absco2 &
    1382    97320960 :                  + wx(2,lay) * cfc11adj(ig) &
    1383   486604800 :                  + wx(3,lay) * cfc12(ig)
    1384   109486080 :             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    18800640 :       do lay = laytrop+1, nlayers
    1391             : 
    1392   164782080 :          do ig = 1, ng6
    1393   291962880 :             taug(lay,ngs5+ig) = 0.0_r8 &
    1394   291962880 :                  + wx(2,lay) * cfc11adj(ig) &
    1395   437944320 :                  + wx(3,lay) * cfc12(ig)
    1396   164229120 :             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    12718080 :       do lay = 1, laytrop
    1455             : 
    1456    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oo3(lay)*colo3(lay)
    1457    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1458    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1459    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    1460    12165120 :          js(lay) = 1 + int(specmult(lay))
    1461    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1462             : 
    1463    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oo3_1(lay)*colo3(lay)
    1464    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1465    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1466    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1467    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    1468    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1469             : 
    1470    12165120 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*colo3(lay)
    1471    12165120 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    1472    12165120 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    1473    12165120 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    1474             : 
    1475    12165120 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    1476    12165120 :          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    12165120 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1482    12165120 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1483    12165120 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    1488             :          endif
    1489             : 
    1490    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colo3(lay)
    1491    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1492    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1493    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1494    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1495    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1496             : 
    1497    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(7) + js(lay)
    1498    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(7) + js1(lay)
    1499             : 
    1500    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
    1501     2489629 :             p = fs(lay) - 1
    1502     2489629 :             p4 = p**4
    1503     2489629 :             fk0 = p4
    1504     2489629 :             fk1 = 1 - p - 2.0_r8*p4
    1505     2489629 :             fk2 = p + p4
    1506     2489629 :             fac000(lay) = fk0*fac00(lay)
    1507     2489629 :             fac100(lay) = fk1*fac00(lay)
    1508     2489629 :             fac200(lay) = fk2*fac00(lay)
    1509     2489629 :             fac010(lay) = fk0*fac10(lay)
    1510     2489629 :             fac110(lay) = fk1*fac10(lay)
    1511     2489629 :             fac210(lay) = fk2*fac10(lay)
    1512     9675491 :          else if (specparm(lay) .gt. 0.875_r8) then
    1513      874144 :             p = -fs(lay) 
    1514      874144 :             p4 = p**4
    1515      874144 :             fk0 = p4
    1516      874144 :             fk1 = 1 - p - 2.0_r8*p4
    1517      874144 :             fk2 = p + p4
    1518      874144 :             fac000(lay) = fk0*fac00(lay)
    1519      874144 :             fac100(lay) = fk1*fac00(lay)
    1520      874144 :             fac200(lay) = fk2*fac00(lay)
    1521      874144 :             fac010(lay) = fk0*fac10(lay)
    1522      874144 :             fac110(lay) = fk1*fac10(lay)
    1523      874144 :             fac210(lay) = fk2*fac10(lay)
    1524             :          else
    1525     8801347 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1526     8801347 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1527     8801347 :             fac100(lay) = fs(lay) * fac00(lay)
    1528     8801347 :             fac110(lay) = fs(lay) * fac10(lay)
    1529             :          endif
    1530    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
    1531     2489629 :             p = fs1(lay) - 1
    1532     2489629 :             p4 = p**4
    1533     2489629 :             fk0 = p4
    1534     2489629 :             fk1 = 1 - p - 2.0_r8*p4
    1535     2489629 :             fk2 = p + p4
    1536     2489629 :             fac001(lay) = fk0*fac01(lay)
    1537     2489629 :             fac101(lay) = fk1*fac01(lay)
    1538     2489629 :             fac201(lay) = fk2*fac01(lay)
    1539     2489629 :             fac011(lay) = fk0*fac11(lay)
    1540     2489629 :             fac111(lay) = fk1*fac11(lay)
    1541     2489629 :             fac211(lay) = fk2*fac11(lay)
    1542     9675491 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1543     1946570 :             p = -fs1(lay) 
    1544     1946570 :             p4 = p**4
    1545     1946570 :             fk0 = p4
    1546     1946570 :             fk1 = 1 - p - 2.0_r8*p4
    1547     1946570 :             fk2 = p + p4
    1548     1946570 :             fac001(lay) = fk0*fac01(lay)
    1549     1946570 :             fac101(lay) = fk1*fac01(lay)
    1550     1946570 :             fac201(lay) = fk2*fac01(lay)
    1551     1946570 :             fac011(lay) = fk0*fac11(lay)
    1552     1946570 :             fac111(lay) = fk1*fac11(lay)
    1553     1946570 :             fac211(lay) = fk2*fac11(lay)
    1554             :          else
    1555     7728921 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1556     7728921 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1557     7728921 :             fac101(lay) = fs1(lay) * fac01(lay)
    1558     7728921 :             fac111(lay) = fs1(lay) * fac11(lay)
    1559             :          endif
    1560             : 
    1561   158699520 :          do ig = 1, ng7
    1562   437944320 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1563   437944320 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1564   291962880 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1565   291962880 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1566   437944320 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    1567   437944320 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    1568   145981440 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    1569   145981440 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    1570   145981440 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    1571             : 
    1572   145981440 :             if (specparm(lay) .lt. 0.125_r8) then
    1573             :                tau_major = speccomb(lay) * &
    1574    29875548 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1575    29875548 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1576    29875548 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1577    29875548 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1578    29875548 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1579   179253288 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1580   116105892 :             else if (specparm(lay) .gt. 0.875_r8) then
    1581             :                tau_major = speccomb(lay) * &
    1582    10489728 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1583    10489728 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1584    10489728 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1585    10489728 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1586    10489728 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1587    62938368 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1588             :             else
    1589             :                tau_major = speccomb(lay) * &
    1590   105616164 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1591   105616164 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1592   105616164 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1593   422464656 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1594             :             endif
    1595             : 
    1596   145981440 :             if (specparm1(lay) .lt. 0.125_r8) then
    1597             :                tau_major1 = speccomb1(lay) * &
    1598    11263824 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    1599    11263824 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1600    11263824 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1601    11263824 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1602    11263824 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1603    67582944 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1604   134717616 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1605             :                tau_major1 = speccomb1(lay) * &
    1606    23358840 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1607    23358840 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1608    23358840 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1609    23358840 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1610    23358840 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    1611   140153040 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    1612             :             else
    1613             :                tau_major1 = speccomb1(lay) * &
    1614   111358776 :                     (fac001(lay) * absa(ind1(lay),ig) +  &
    1615   111358776 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1616   111358776 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1617   445435104 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    1618             :             endif
    1619             : 
    1620   145981440 :             taug(lay,ngs6+ig) = tau_major + tau_major1 &
    1621             :                  + tauself + taufor &
    1622   145981440 :                  + adjcolco2(lay)*absco2
    1623   291962880 :             fracs(lay,ngs6+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    1624   304128000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    1625             :          enddo
    1626             :       enddo
    1627             : 
    1628             : ! Upper atmosphere loop
    1629    18800640 :       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    18247680 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1635    18247680 :          ratco2(lay) = 1.e20*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1636    18247680 :          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    18247680 :             adjcolco2(lay) = colco2(lay)
    1641             :          endif
    1642             : 
    1643    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(7) + 1
    1644    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(7) + 1
    1645             : 
    1646   237219840 :          do ig = 1, ng7
    1647   437944320 :             absco2 = kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1648   656916480 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig))
    1649   218972160 :             taug(lay,ngs6+ig) = colo3(lay) * &
    1650   437944320 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1651   437944320 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1652   437944320 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1653   437944320 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1654  1094860800 :                  + adjcolco2(lay) * absco2
    1655   237219840 :             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    18247680 :          taug(lay,ngs6+6)=taug(lay,ngs6+6)*0.92_r8
    1662    18247680 :          taug(lay,ngs6+7)=taug(lay,ngs6+7)*0.88_r8
    1663    18247680 :          taug(lay,ngs6+8)=taug(lay,ngs6+8)*1.07_r8
    1664    18247680 :          taug(lay,ngs6+9)=taug(lay,ngs6+9)*1.1_r8
    1665    18247680 :          taug(lay,ngs6+10)=taug(lay,ngs6+10)*0.99_r8
    1666    18800640 :          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    12718080 :       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    12165120 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    1716    12165120 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1717    12165120 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    1722             :          endif
    1723             : 
    1724    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(8) + 1
    1725    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(8) + 1
    1726             : 
    1727   110039040 :          do ig = 1, ng8
    1728   291962880 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    1729   291962880 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1730   194641920 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1731   194641920 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    1732   194641920 :             absco2 =  (ka_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1733   194641920 :                  (ka_mco2(indminor(lay)+1,ig) - ka_mco2(indminor(lay),ig)))
    1734             :             abso3 =  (ka_mo3(indminor(lay),ig) + minorfrac(lay) * &
    1735    97320960 :                  (ka_mo3(indminor(lay)+1,ig) - ka_mo3(indminor(lay),ig)))
    1736             :             absn2o =  (ka_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1737    97320960 :                  (ka_mn2o(indminor(lay)+1,ig) - ka_mn2o(indminor(lay),ig)))
    1738    97320960 :             taug(lay,ngs7+ig) = colh2o(lay) * &
    1739   194641920 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    1740   194641920 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    1741   194641920 :                  fac01(lay) * absa(ind1(lay),ig) +  &
    1742   194641920 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    1743             :                  + tauself + taufor &
    1744             :                  + adjcolco2(lay)*absco2 &
    1745    97320960 :                  + colo3(lay) * abso3 &
    1746    97320960 :                  + coln2o(lay) * absn2o &
    1747    97320960 :                  + wx(3,lay) * cfc12(ig) &
    1748   486604800 :                  + wx(4,lay) * cfc22adj(ig)
    1749   109486080 :             fracs(lay,ngs7+ig) = fracrefa(ig)
    1750             :          enddo
    1751             :       enddo
    1752             : 
    1753             : ! Upper atmosphere loop
    1754    18800640 :       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    18247680 :          chi_co2(lay) = colco2(lay)/coldry(lay)
    1760    18247680 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/chi_mls(2,jp(lay)+1)
    1761    18247680 :          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    18247680 :             adjcolco2(lay) = colco2(lay)
    1766             :          endif
    1767             : 
    1768    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(8) + 1
    1769    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(8) + 1
    1770             : 
    1771   164782080 :          do ig = 1, ng8
    1772   291962880 :             absco2 =  (kb_mco2(indminor(lay),ig) + minorfrac(lay) * &
    1773   437944320 :                  (kb_mco2(indminor(lay)+1,ig) - kb_mco2(indminor(lay),ig)))
    1774             :             absn2o =  (kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    1775   145981440 :                  (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig)))
    1776   145981440 :             taug(lay,ngs7+ig) = colo3(lay) * &
    1777   291962880 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    1778   291962880 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    1779   291962880 :                  fac01(lay) * absb(ind1(lay),ig) + &
    1780   291962880 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    1781             :                  + adjcolco2(lay)*absco2 &
    1782   145981440 :                  + coln2o(lay)*absn2o & 
    1783   145981440 :                  + wx(3,lay) * cfc12(ig) &
    1784   729907200 :                  + wx(4,lay) * cfc22adj(ig)
    1785   164229120 :             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    12718080 :       do lay = 1, laytrop
    1844             : 
    1845    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    1846    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    1847    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    1848    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    1849    12165120 :          js(lay) = 1 + int(specmult(lay))
    1850    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    1851             : 
    1852    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    1853    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    1854    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    1855    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    1856    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    1857    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    1858             : 
    1859    12165120 :          speccomb_mn2o(lay) = colh2o(lay) + refrat_m_a*colch4(lay)
    1860    12165120 :          specparm_mn2o(lay) = colh2o(lay)/speccomb_mn2o(lay)
    1861    12165120 :          if (specparm_mn2o(lay) .ge. oneminus) specparm_mn2o(lay) = oneminus
    1862    12165120 :          specmult_mn2o(lay) = 8._r8*specparm_mn2o(lay)
    1863    12165120 :          jmn2o(lay) = 1 + int(specmult_mn2o(lay))
    1864    12165120 :          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    12165120 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    1870    12165120 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    1871    12165120 :          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    12165120 :             adjcoln2o(lay) = coln2o(lay)
    1876             :          endif
    1877             : 
    1878    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    1879    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    1880    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    1881    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    1882    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    1883    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    1884             : 
    1885    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(9) + js(lay)
    1886    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(9) + js1(lay)
    1887             : 
    1888    12165120 :          if (specparm(lay) .lt. 0.125_r8) then
    1889     2156503 :             p = fs(lay) - 1
    1890     2156503 :             p4 = p**4
    1891     2156503 :             fk0 = p4
    1892     2156503 :             fk1 = 1 - p - 2.0_r8*p4
    1893     2156503 :             fk2 = p + p4
    1894     2156503 :             fac000(lay) = fk0*fac00(lay)
    1895     2156503 :             fac100(lay) = fk1*fac00(lay)
    1896     2156503 :             fac200(lay) = fk2*fac00(lay)
    1897     2156503 :             fac010(lay) = fk0*fac10(lay)
    1898     2156503 :             fac110(lay) = fk1*fac10(lay)
    1899     2156503 :             fac210(lay) = fk2*fac10(lay)
    1900    10008617 :          else if (specparm(lay) .gt. 0.875_r8) then
    1901       29999 :             p = -fs(lay) 
    1902       29999 :             p4 = p**4
    1903       29999 :             fk0 = p4
    1904       29999 :             fk1 = 1 - p - 2.0_r8*p4
    1905       29999 :             fk2 = p + p4
    1906       29999 :             fac000(lay) = fk0*fac00(lay)
    1907       29999 :             fac100(lay) = fk1*fac00(lay)
    1908       29999 :             fac200(lay) = fk2*fac00(lay)
    1909       29999 :             fac010(lay) = fk0*fac10(lay)
    1910       29999 :             fac110(lay) = fk1*fac10(lay)
    1911       29999 :             fac210(lay) = fk2*fac10(lay)
    1912             :          else
    1913     9978618 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    1914     9978618 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    1915     9978618 :             fac100(lay) = fs(lay) * fac00(lay)
    1916     9978618 :             fac110(lay) = fs(lay) * fac10(lay)
    1917             :          endif
    1918             : 
    1919    12165120 :          if (specparm1(lay) .lt. 0.125_r8) then
    1920      845664 :             p = fs1(lay) - 1
    1921      845664 :             p4 = p**4
    1922      845664 :             fk0 = p4
    1923      845664 :             fk1 = 1 - p - 2.0_r8*p4
    1924      845664 :             fk2 = p + p4
    1925      845664 :             fac001(lay) = fk0*fac01(lay)
    1926      845664 :             fac101(lay) = fk1*fac01(lay)
    1927      845664 :             fac201(lay) = fk2*fac01(lay)
    1928      845664 :             fac011(lay) = fk0*fac11(lay)
    1929      845664 :             fac111(lay) = fk1*fac11(lay)
    1930      845664 :             fac211(lay) = fk2*fac11(lay)
    1931    11319456 :          else if (specparm1(lay) .gt. 0.875_r8) then
    1932      365328 :             p = -fs1(lay) 
    1933      365328 :             p4 = p**4
    1934      365328 :             fk0 = p4
    1935      365328 :             fk1 = 1 - p - 2.0_r8*p4
    1936      365328 :             fk2 = p + p4
    1937      365328 :             fac001(lay) = fk0*fac01(lay)
    1938      365328 :             fac101(lay) = fk1*fac01(lay)
    1939      365328 :             fac201(lay) = fk2*fac01(lay)
    1940      365328 :             fac011(lay) = fk0*fac11(lay)
    1941      365328 :             fac111(lay) = fk1*fac11(lay)
    1942      365328 :             fac211(lay) = fk2*fac11(lay)
    1943             :          else
    1944    10954128 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    1945    10954128 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    1946    10954128 :             fac101(lay) = fs1(lay) * fac01(lay)
    1947    10954128 :             fac111(lay) = fs1(lay) * fac11(lay)
    1948             :          endif
    1949             : 
    1950   158699520 :          do ig = 1, ng9
    1951   437944320 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    1952   437944320 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    1953   291962880 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    1954   291962880 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    1955   437944320 :             n2om1 = ka_mn2o(jmn2o(lay),indminor(lay),ig) + fmn2o(lay) * &
    1956   437944320 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay),ig) - ka_mn2o(jmn2o(lay),indminor(lay),ig))
    1957   145981440 :             n2om2 = ka_mn2o(jmn2o(lay),indminor(lay)+1,ig) + fmn2o(lay) * &
    1958   145981440 :                  (ka_mn2o(jmn2o(lay)+1,indminor(lay)+1,ig) - ka_mn2o(jmn2o(lay),indminor(lay)+1,ig))
    1959   145981440 :             absn2o = n2om1 + minorfrac(lay) * (n2om2 - n2om1)
    1960             : 
    1961   145981440 :             if (specparm(lay) .lt. 0.125_r8) then
    1962             :                tau_major = speccomb(lay) * &
    1963    25878036 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1964    25878036 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1965    25878036 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    1966    25878036 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1967    25878036 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    1968   155268216 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    1969   120103404 :             else if (specparm(lay) .gt. 0.875_r8) then
    1970             :                tau_major = speccomb(lay) * &
    1971      359988 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    1972      359988 :                     fac100(lay) * absa(ind0(lay),ig) + &
    1973      359988 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    1974      359988 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    1975      359988 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    1976     2159928 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    1977             :             else
    1978             :                tau_major = speccomb(lay) * &
    1979   119743416 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    1980   119743416 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    1981   119743416 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    1982   478973664 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    1983             :             endif
    1984             : 
    1985   145981440 :             if (specparm1(lay) .lt. 0.125_r8) then
    1986             :                tau_major1 = speccomb1(lay) * &
    1987    10147968 :                     (fac001(lay) * absa(ind1(lay),ig) + & 
    1988    10147968 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    1989    10147968 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    1990    10147968 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    1991    10147968 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    1992    60887808 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    1993   135833472 :             else if (specparm1(lay) .gt. 0.875_r8) then
    1994             :                tau_major1 = speccomb1(lay) * &
    1995     4383936 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    1996     4383936 :                     fac101(lay) * absa(ind1(lay),ig) + &
    1997     4383936 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    1998     4383936 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    1999     4383936 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2000    26303616 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2001             :             else
    2002             :                tau_major1 = speccomb1(lay) * &
    2003   131449536 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2004   131449536 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2005   131449536 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2006   525798144 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2007             :             endif
    2008             : 
    2009   145981440 :             taug(lay,ngs8+ig) = tau_major + tau_major1 &
    2010             :                  + tauself + taufor &
    2011   145981440 :                  + adjcoln2o(lay)*absn2o
    2012   291962880 :             fracs(lay,ngs8+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2013   304128000 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2014             :          enddo
    2015             :       enddo
    2016             : 
    2017             : ! Upper atmosphere loop
    2018    18800640 :       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    18247680 :          chi_n2o(lay) = coln2o(lay)/(coldry(lay))
    2024    18247680 :          ratn2o(lay) = 1.e20_r8*chi_n2o(lay)/chi_mls(4,jp(lay)+1)
    2025    18247680 :          if (ratn2o(lay) .gt. 1.5_r8) then
    2026     4897891 :             adjfac(lay) = 0.5_r8+(ratn2o(lay)-0.5_r8)**0.65_r8
    2027     4897891 :             adjcoln2o(lay) = adjfac(lay)*chi_mls(4,jp(lay)+1)*coldry(lay)*1.e-20_r8
    2028             :          else
    2029    13349789 :             adjcoln2o(lay) = coln2o(lay)
    2030             :          endif
    2031             : 
    2032    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(9) + 1
    2033    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(9) + 1
    2034             : 
    2035   237772800 :          do ig = 1, ng9
    2036   437944320 :             absn2o = kb_mn2o(indminor(lay),ig) + minorfrac(lay) * &
    2037   656916480 :                 (kb_mn2o(indminor(lay)+1,ig) - kb_mn2o(indminor(lay),ig))
    2038   218972160 :             taug(lay,ngs8+ig) = colch4(lay) * &
    2039   437944320 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2040   437944320 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2041   437944320 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2042   437944320 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2043  1094860800 :                  + adjcoln2o(lay)*absn2o
    2044   237219840 :             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    12718080 :       do lay = 1, laytrop
    2076    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(10) + 1
    2077    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(10) + 1
    2078             : 
    2079    85708800 :          do ig = 1, ng10
    2080   218972160 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2081   218972160 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2082   145981440 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2083   145981440 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2084    72990720 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2085   145981440 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2086   145981440 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2087   145981440 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2088   145981440 :                  fac11(lay) * absa(ind1(lay)+1,ig))  &
    2089   364953600 :                  + tauself + taufor
    2090    85155840 :             fracs(lay,ngs9+ig) = fracrefa(ig)
    2091             :          enddo
    2092             :       enddo
    2093             : 
    2094             : ! Upper atmosphere loop
    2095    18800640 :       do lay = laytrop+1, nlayers
    2096    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(10) + 1
    2097    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(10) + 1
    2098             : 
    2099   128286720 :          do ig = 1, ng10
    2100   328458240 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2101   328458240 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2102   109486080 :             taug(lay,ngs9+ig) = colh2o(lay) * &
    2103   218972160 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2104   218972160 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2105   218972160 :                  fac01(lay) * absb(ind1(lay),ig) +  &
    2106   218972160 :                  fac11(lay) * absb(ind1(lay)+1,ig)) &
    2107   547430400 :                  + taufor
    2108   127733760 :             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    12718080 :       do lay = 1, laytrop
    2145    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(11) + 1
    2146    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(11) + 1
    2147    12165120 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2148   110039040 :          do ig = 1, ng11
    2149   291962880 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2150   291962880 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2151   194641920 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2152   194641920 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig)))
    2153   194641920 :             tauo2 =  scaleo2(lay) * (ka_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2154   194641920 :                  (ka_mo2(indminor(lay)+1,ig) - ka_mo2(indminor(lay),ig)))
    2155    97320960 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2156   194641920 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2157   194641920 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2158   194641920 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2159   194641920 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2160             :                  + tauself + taufor &
    2161   486604800 :                  + tauo2
    2162   109486080 :             fracs(lay,ngs10+ig) = fracrefa(ig)
    2163             :          enddo
    2164             :       enddo
    2165             : 
    2166             : ! Upper atmosphere loop
    2167    18800640 :       do lay = laytrop+1, nlayers
    2168    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(11) + 1
    2169    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(11) + 1
    2170    18247680 :          scaleo2(lay) = colo2(lay)*scaleminor(lay)
    2171   164782080 :          do ig = 1, ng11
    2172   437944320 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2173   437944320 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2174   291962880 :             tauo2 =  scaleo2(lay) * (kb_mo2(indminor(lay),ig) + minorfrac(lay) * &
    2175   291962880 :                  (kb_mo2(indminor(lay)+1,ig) - kb_mo2(indminor(lay),ig)))
    2176   145981440 :             taug(lay,ngs10+ig) = colh2o(lay) * &
    2177   291962880 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2178   291962880 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2179   291962880 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2180   291962880 :                  fac11(lay) * absb(ind1(lay)+1,ig))  &
    2181             :                  + taufor &
    2182   729907200 :                  + tauo2
    2183   164229120 :             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    12718080 :       do lay = 1, laytrop
    2232             : 
    2233    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2oco2(lay)*colco2(lay)
    2234    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2235    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2236    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    2237    12165120 :          js(lay) = 1 + int(specmult(lay))
    2238    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2239             : 
    2240    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2oco2_1(lay)*colco2(lay)
    2241    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2242    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2243    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2244    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    2245    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2246             : 
    2247    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colco2(lay)
    2248    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2249    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2250    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2251    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2252    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2253             : 
    2254    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(12) + js(lay)
    2255    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(12) + js1(lay)
    2256             :       enddo
    2257    12718080 :       do lay = 1, laytrop
    2258    12718080 :          if (specparm(lay) .lt. 0.125_r8) then
    2259     2698634 :             p = fs(lay) - 1
    2260     2698634 :             p4 = p**4
    2261     2698634 :             fk0 = p4
    2262     2698634 :             fk1 = 1 - p - 2.0_r8*p4
    2263     2698634 :             fk2 = p + p4
    2264     2698634 :             fac000(lay) = fk0*fac00(lay)
    2265     2698634 :             fac100(lay) = fk1*fac00(lay)
    2266     2698634 :             fac200(lay) = fk2*fac00(lay)
    2267     2698634 :             fac010(lay) = fk0*fac10(lay)
    2268     2698634 :             fac110(lay) = fk1*fac10(lay)
    2269     2698634 :             fac210(lay) = fk2*fac10(lay)
    2270     9466486 :          else if (specparm(lay) .gt. 0.875_r8) then
    2271       13654 :             p = -fs(lay) 
    2272       13654 :             p4 = p**4
    2273       13654 :             fk0 = p4
    2274       13654 :             fk1 = 1 - p - 2.0_r8*p4
    2275       13654 :             fk2 = p + p4
    2276       13654 :             fac000(lay) = fk0*fac00(lay)
    2277       13654 :             fac100(lay) = fk1*fac00(lay)
    2278       13654 :             fac200(lay) = fk2*fac00(lay)
    2279       13654 :             fac010(lay) = fk0*fac10(lay)
    2280       13654 :             fac110(lay) = fk1*fac10(lay)
    2281       13654 :             fac210(lay) = fk2*fac10(lay)
    2282             :          else
    2283     9452832 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2284     9452832 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2285     9452832 :             fac100(lay) = fs(lay) * fac00(lay)
    2286     9452832 :             fac110(lay) = fs(lay) * fac10(lay)
    2287             :          endif
    2288             : 
    2289             :       enddo
    2290    12718080 :       do lay = 1, laytrop
    2291    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
    2292     1146527 :             p = fs1(lay) - 1
    2293     1146527 :             p4 = p**4
    2294     1146527 :             fk0 = p4
    2295     1146527 :             fk1 = 1 - p - 2.0_r8*p4
    2296     1146527 :             fk2 = p + p4
    2297     1146527 :             fac001(lay) = fk0*fac01(lay)
    2298     1146527 :             fac101(lay) = fk1*fac01(lay)
    2299     1146527 :             fac201(lay) = fk2*fac01(lay)
    2300     1146527 :             fac011(lay) = fk0*fac11(lay)
    2301     1146527 :             fac111(lay) = fk1*fac11(lay)
    2302     1146527 :             fac211(lay) = fk2*fac11(lay)
    2303    11018593 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2304      255010 :             p = -fs1(lay) 
    2305      255010 :             p4 = p**4
    2306      255010 :             fk0 = p4
    2307      255010 :             fk1 = 1 - p - 2.0_r8*p4
    2308      255010 :             fk2 = p + p4
    2309      255010 :             fac001(lay) = fk0*fac01(lay)
    2310      255010 :             fac101(lay) = fk1*fac01(lay)
    2311      255010 :             fac201(lay) = fk2*fac01(lay)
    2312      255010 :             fac011(lay) = fk0*fac11(lay)
    2313      255010 :             fac111(lay) = fk1*fac11(lay)
    2314      255010 :             fac211(lay) = fk2*fac11(lay)
    2315             :          else
    2316    10763583 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2317    10763583 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2318    10763583 :             fac101(lay) = fs1(lay) * fac01(lay)
    2319    10763583 :             fac111(lay) = fs1(lay) * fac11(lay)
    2320             :          endif
    2321             :       enddo
    2322    12718080 :       do lay = 1, laytrop
    2323             : 
    2324   110039040 :          do ig = 1, ng12
    2325   291962880 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2326   389283840 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2327   194641920 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2328   194641920 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2329             : 
    2330    97320960 :             if (specparm(lay) .lt. 0.125_r8) then
    2331             :                tau_major = speccomb(lay) * &
    2332    21589072 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2333    21589072 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2334    21589072 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2335    21589072 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2336    21589072 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2337   129534432 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2338    75731888 :             else if (specparm(lay) .gt. 0.875_r8) then
    2339             :                tau_major = speccomb(lay) * &
    2340      109232 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2341      109232 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2342      109232 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2343      109232 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2344      109232 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2345      655392 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2346             :             else
    2347             :                tau_major = speccomb(lay) * &
    2348    75622656 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2349    75622656 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2350    75622656 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2351   302490624 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2352             :             endif
    2353             : 
    2354    97320960 :             if (specparm1(lay) .lt. 0.125_r8) then
    2355             :                tau_major1 = speccomb1(lay) * &
    2356     9172216 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2357     9172216 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2358     9172216 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2359     9172216 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2360     9172216 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2361    55033296 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2362    88148744 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2363             :                tau_major1 = speccomb1(lay) * &
    2364     2040080 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2365     2040080 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2366     2040080 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2367     2040080 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2368     2040080 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2369    12240480 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2370             :             else
    2371             :                tau_major1 = speccomb1(lay) * &
    2372    86108664 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2373    86108664 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2374    86108664 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2375   344434656 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2376             :             endif
    2377             : 
    2378    97320960 :             taug(lay,ngs11+ig) = tau_major + tau_major1 &
    2379    97320960 :                  + tauself + taufor
    2380   194641920 :             fracs(lay,ngs11+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2381   206807040 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2382             :          enddo
    2383             :       enddo
    2384             : 
    2385             : ! Upper atmosphere loop
    2386    18800640 :       do lay = laytrop+1, nlayers
    2387   164782080 :          do ig = 1, ng12
    2388   145981440 :             taug(lay,ngs11+ig) = 0.0_r8
    2389   164229120 :             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    12718080 :       do lay = 1, laytrop
    2453             : 
    2454    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2on2o(lay)*coln2o(lay)
    2455    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2456    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2457    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    2458    12165120 :          js(lay) = 1 + int(specmult(lay))
    2459    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2460             : 
    2461    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2on2o_1(lay)*coln2o(lay)
    2462    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    2463    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2464    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2465    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    2466    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2467             : 
    2468    12165120 :          speccomb_mco2(lay) = colh2o(lay) + refrat_m_a*coln2o(lay)
    2469    12165120 :          specparm_mco2(lay) = colh2o(lay)/speccomb_mco2(lay)
    2470    12165120 :          if (specparm_mco2(lay) .ge. oneminus) specparm_mco2(lay) = oneminus
    2471    12165120 :          specmult_mco2(lay) = 8._r8*specparm_mco2(lay)
    2472    12165120 :          jmco2(lay) = 1 + int(specmult_mco2(lay))
    2473    12165120 :          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    12165120 :          chi_co2(lay) = colco2(lay)/(coldry(lay))
    2479    12718080 :          ratco2(lay) = 1.e20_r8*chi_co2(lay)/3.55e-4_r8
    2480             :       enddo
    2481    12718080 :       do lay = 1, laytrop
    2482    12718080 :          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    12165120 :             adjcolco2(lay) = colco2(lay)
    2487             :          endif
    2488             : 
    2489             :       enddo
    2490    12718080 :       do lay = 1, laytrop
    2491    12165120 :          speccomb_mco(lay) = colh2o(lay) + refrat_m_a3*coln2o(lay)
    2492    12165120 :          specparm_mco(lay) = colh2o(lay)/speccomb_mco(lay)
    2493    12165120 :          if (specparm_mco(lay) .ge. oneminus) specparm_mco(lay) = oneminus
    2494    12165120 :          specmult_mco(lay) = 8._r8*specparm_mco(lay)
    2495    12165120 :          jmco(lay) = 1 + int(specmult_mco(lay))
    2496    12165120 :          fmco(lay) = mod(specmult_mco(lay),1.0_r8)
    2497             : 
    2498    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*coln2o(lay)
    2499    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    2500    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2501    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2502    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2503    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2504             : 
    2505    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(13) + js(lay)
    2506    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(13) + js1(lay)
    2507             :       enddo
    2508    12718080 :       do lay = 1, laytrop
    2509             : 
    2510    12718080 :          if (specparm(lay) .lt. 0.125_r8) then
    2511     2785653 :             p = fs(lay) - 1
    2512     2785653 :             p4 = p**4
    2513     2785653 :             fk0 = p4
    2514     2785653 :             fk1 = 1 - p - 2.0_r8*p4
    2515     2785653 :             fk2 = p + p4
    2516     2785653 :             fac000(lay) = fk0*fac00(lay)
    2517     2785653 :             fac100(lay) = fk1*fac00(lay)
    2518     2785653 :             fac200(lay) = fk2*fac00(lay)
    2519     2785653 :             fac010(lay) = fk0*fac10(lay)
    2520     2785653 :             fac110(lay) = fk1*fac10(lay)
    2521     2785653 :             fac210(lay) = fk2*fac10(lay)
    2522     9379467 :          else if (specparm(lay) .gt. 0.875_r8) then
    2523        5828 :             p = -fs(lay) 
    2524        5828 :             p4 = p**4
    2525        5828 :             fk0 = p4
    2526        5828 :             fk1 = 1 - p - 2.0_r8*p4
    2527        5828 :             fk2 = p + p4
    2528        5828 :             fac000(lay) = fk0*fac00(lay)
    2529        5828 :             fac100(lay) = fk1*fac00(lay)
    2530        5828 :             fac200(lay) = fk2*fac00(lay)
    2531        5828 :             fac010(lay) = fk0*fac10(lay)
    2532        5828 :             fac110(lay) = fk1*fac10(lay)
    2533        5828 :             fac210(lay) = fk2*fac10(lay)
    2534             :          else
    2535     9373639 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2536     9373639 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2537     9373639 :             fac100(lay) = fs(lay) * fac00(lay)
    2538     9373639 :             fac110(lay) = fs(lay) * fac10(lay)
    2539             :          endif
    2540             : 
    2541             :       enddo
    2542    12718080 :       do lay = 1, laytrop
    2543    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
    2544     1206242 :             p = fs1(lay) - 1
    2545     1206242 :             p4 = p**4
    2546     1206242 :             fk0 = p4
    2547     1206242 :             fk1 = 1 - p - 2.0_r8*p4
    2548     1206242 :             fk2 = p + p4
    2549     1206242 :             fac001(lay) = fk0*fac01(lay)
    2550     1206242 :             fac101(lay) = fk1*fac01(lay)
    2551     1206242 :             fac201(lay) = fk2*fac01(lay)
    2552     1206242 :             fac011(lay) = fk0*fac11(lay)
    2553     1206242 :             fac111(lay) = fk1*fac11(lay)
    2554     1206242 :             fac211(lay) = fk2*fac11(lay)
    2555    10958878 :          else if (specparm1(lay) .gt. 0.875_r8) then
    2556      209939 :             p = -fs1(lay) 
    2557      209939 :             p4 = p**4
    2558      209939 :             fk0 = p4
    2559      209939 :             fk1 = 1 - p - 2.0_r8*p4
    2560      209939 :             fk2 = p + p4
    2561      209939 :             fac001(lay) = fk0*fac01(lay)
    2562      209939 :             fac101(lay) = fk1*fac01(lay)
    2563      209939 :             fac201(lay) = fk2*fac01(lay)
    2564      209939 :             fac011(lay) = fk0*fac11(lay)
    2565      209939 :             fac111(lay) = fk1*fac11(lay)
    2566      209939 :             fac211(lay) = fk2*fac11(lay)
    2567             :          else
    2568    10748939 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2569    10748939 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2570    10748939 :             fac101(lay) = fs1(lay) * fac01(lay)
    2571    10748939 :             fac111(lay) = fs1(lay) * fac11(lay)
    2572             :          endif
    2573             :       enddo
    2574    12718080 :       do lay = 1, laytrop
    2575             : 
    2576    61378560 :          do ig = 1, ng13
    2577   145981440 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2578   194641920 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2579    97320960 :             taufor = forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2580    97320960 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2581   145981440 :             co2m1 = ka_mco2(jmco2(lay),indminor(lay),ig) + fmco2(lay) * &
    2582   145981440 :                  (ka_mco2(jmco2(lay)+1,indminor(lay),ig) - ka_mco2(jmco2(lay),indminor(lay),ig))
    2583    48660480 :             co2m2 = ka_mco2(jmco2(lay),indminor(lay)+1,ig) + fmco2(lay) * &
    2584    48660480 :                  (ka_mco2(jmco2(lay)+1,indminor(lay)+1,ig) - ka_mco2(jmco2(lay),indminor(lay)+1,ig))
    2585    48660480 :             absco2 = co2m1 + minorfrac(lay) * (co2m2 - co2m1)
    2586    48660480 :             com1 = ka_mco(jmco(lay),indminor(lay),ig) + fmco(lay) * &
    2587    97320960 :                  (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    48660480 :                  (ka_mco(jmco(lay)+1,indminor(lay)+1,ig) - ka_mco(jmco(lay),indminor(lay)+1,ig))
    2590    48660480 :             absco = com1 + minorfrac(lay) * (com2 - com1)
    2591             : 
    2592    48660480 :             if (specparm(lay) .lt. 0.125_r8) then
    2593             :                tau_major = speccomb(lay) * &
    2594    11142612 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2595    11142612 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2596    11142612 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    2597    11142612 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2598    11142612 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    2599    66855672 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    2600    37517868 :             else if (specparm(lay) .gt. 0.875_r8) then
    2601             :                tau_major = speccomb(lay) * &
    2602       23312 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    2603       23312 :                     fac100(lay) * absa(ind0(lay),ig) + &
    2604       23312 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    2605       23312 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    2606       23312 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    2607      139872 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    2608             :             else
    2609             :                tau_major = speccomb(lay) * &
    2610    37494556 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2611    37494556 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2612    37494556 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2613   149978224 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2614             :             endif
    2615             : 
    2616    48660480 :             if (specparm1(lay) .lt. 0.125_r8) then
    2617             :                tau_major1 = speccomb1(lay) * &
    2618     4824968 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2619     4824968 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2620     4824968 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    2621     4824968 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2622     4824968 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    2623    28949808 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    2624    43835512 :             else if (specparm1(lay) .gt. 0.875_r8) then
    2625             :                tau_major1 = speccomb1(lay) * &
    2626      839756 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    2627      839756 :                     fac101(lay) * absa(ind1(lay),ig) + &
    2628      839756 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    2629      839756 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    2630      839756 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    2631     5038536 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    2632             :             else
    2633             :                tau_major1 = speccomb1(lay) * &
    2634    42995756 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2635    42995756 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2636    42995756 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2637   171983024 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2638             :             endif
    2639             : 
    2640    48660480 :             taug(lay,ngs12+ig) = tau_major + tau_major1 &
    2641             :                  + tauself + taufor &
    2642             :                  + adjcolco2(lay)*absco2 &
    2643    97320960 :                  + colco(lay)*absco
    2644    97320960 :             fracs(lay,ngs12+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2645   109486080 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2646             :          enddo
    2647             :       enddo
    2648             : 
    2649             : ! Upper atmosphere loop
    2650    18800640 :       do lay = laytrop+1, nlayers
    2651    91791360 :          do ig = 1, ng13
    2652   218972160 :             abso3 = kb_mo3(indminor(lay),ig) + minorfrac(lay) * &
    2653   291962880 :                  (kb_mo3(indminor(lay)+1,ig) - kb_mo3(indminor(lay),ig))
    2654    72990720 :             taug(lay,ngs12+ig) = colo3(lay)*abso3
    2655    91238400 :             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    12718080 :       do lay = 1, laytrop
    2687    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(14) + 1
    2688    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(14) + 1
    2689    37048320 :          do ig = 1, ng14
    2690    72990720 :             tauself = selffac(lay) * (selfref(indself(lay),ig) + selffrac(lay) * &
    2691    72990720 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2692    48660480 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2693    48660480 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2694    24330240 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2695    48660480 :                  (fac00(lay) * absa(ind0(lay),ig) + &
    2696    48660480 :                  fac10(lay) * absa(ind0(lay)+1,ig) + &
    2697    48660480 :                  fac01(lay) * absa(ind1(lay),ig) + &
    2698    48660480 :                  fac11(lay) * absa(ind1(lay)+1,ig)) &
    2699   121651200 :                  + tauself + taufor
    2700    36495360 :             fracs(lay,ngs13+ig) = fracrefa(ig)
    2701             :          enddo
    2702             :       enddo
    2703             : 
    2704             : ! Upper atmosphere loop
    2705    18800640 :       do lay = laytrop+1, nlayers
    2706    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(14) + 1
    2707    18247680 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(14) + 1
    2708    55296000 :          do ig = 1, ng14
    2709    36495360 :             taug(lay,ngs13+ig) = colco2(lay) * &
    2710   109486080 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    2711    72990720 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    2712    72990720 :                  fac01(lay) * absb(ind1(lay),ig) + &
    2713   218972160 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    2714    54743040 :             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    12718080 :       do lay = 1, laytrop
    2770             : 
    2771    12165120 :          speccomb(lay) = coln2o(lay) + rat_n2oco2(lay)*colco2(lay)
    2772    12165120 :          specparm(lay) = coln2o(lay)/speccomb(lay)
    2773    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2774    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    2775    12165120 :          js(lay) = 1 + int(specmult(lay))
    2776    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2777             : 
    2778    12165120 :          speccomb1(lay) = coln2o(lay) + rat_n2oco2_1(lay)*colco2(lay)
    2779    12165120 :          specparm1(lay) = coln2o(lay)/speccomb1(lay)
    2780    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    2781    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    2782    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    2783    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    2784             : 
    2785    12165120 :          speccomb_mn2(lay) = coln2o(lay) + refrat_m_a*colco2(lay)
    2786    12165120 :          specparm_mn2(lay) = coln2o(lay)/speccomb_mn2(lay)
    2787    12165120 :          if (specparm_mn2(lay) .ge. oneminus) specparm_mn2(lay) = oneminus
    2788    12165120 :          specmult_mn2(lay) = 8._r8*specparm_mn2(lay)
    2789    12165120 :          jmn2(lay) = 1 + int(specmult_mn2(lay))
    2790    12165120 :          fmn2(lay) = mod(specmult_mn2(lay),1.0_r8)
    2791             : 
    2792    12165120 :          speccomb_planck(lay) = coln2o(lay)+refrat_planck_a*colco2(lay)
    2793    12165120 :          specparm_planck(lay) = coln2o(lay)/speccomb_planck(lay)
    2794    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    2795    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    2796    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    2797    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    2798             : 
    2799    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(15) + js(lay)
    2800    12165120 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(15) + js1(lay)
    2801             :          
    2802    12718080 :          scalen2(lay) = colbrd(lay)*scaleminor(lay)
    2803             :       enddo
    2804    12718080 :       do lay = 1, laytrop
    2805    12718080 :          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    12165120 :          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    12165120 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    2831    12165120 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    2832    12165120 :             fac100(lay) = fs(lay) * fac00(lay)
    2833    12165120 :             fac110(lay) = fs(lay) * fac10(lay)
    2834             :          endif
    2835             :       enddo
    2836    12718080 :       do lay = 1, laytrop
    2837    12718080 :          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    12165120 :          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    12165120 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    2863    12165120 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    2864    12165120 :             fac101(lay) = fs1(lay) * fac01(lay)
    2865    12165120 :             fac111(lay) = fs1(lay) * fac11(lay)
    2866             :          endif
    2867             : 
    2868             :       enddo
    2869    12718080 :       do lay = 1, laytrop
    2870    37048320 :          do ig = 1, ng15
    2871    72990720 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    2872    97320960 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    2873    48660480 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    2874    48660480 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    2875    72990720 :             n2m1 = ka_mn2(jmn2(lay),indminor(lay),ig) + fmn2(lay) * &
    2876    72990720 :                  (ka_mn2(jmn2(lay)+1,indminor(lay),ig) - ka_mn2(jmn2(lay),indminor(lay),ig))
    2877    24330240 :             n2m2 = ka_mn2(jmn2(lay),indminor(lay)+1,ig) + fmn2(lay) * &
    2878    24330240 :                  (ka_mn2(jmn2(lay)+1,indminor(lay)+1,ig) - ka_mn2(jmn2(lay),indminor(lay)+1,ig))
    2879    24330240 :             taun2 = scalen2(lay) * (n2m1 + minorfrac(lay) * (n2m2 - n2m1))
    2880             : 
    2881    24330240 :             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    24330240 :             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    24330240 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    2900    24330240 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    2901    24330240 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    2902    97320960 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    2903             :             endif 
    2904             : 
    2905    24330240 :             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    24330240 :             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    24330240 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    2924    24330240 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    2925    24330240 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    2926    97320960 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    2927             :             endif
    2928             : 
    2929    24330240 :             taug(lay,ngs14+ig) = tau_major + tau_major1 &
    2930             :                  + tauself + taufor &
    2931    24330240 :                  + taun2
    2932    48660480 :             fracs(lay,ngs14+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    2933    60825600 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    2934             :          enddo
    2935             :       enddo
    2936             : 
    2937             : ! Upper atmosphere loop
    2938    18800640 :       do lay = laytrop+1, nlayers
    2939    55296000 :          do ig = 1, ng15
    2940    36495360 :             taug(lay,ngs14+ig) = 0.0_r8
    2941    54743040 :             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    12718080 :       do lay = 1, laytrop
    2990             : 
    2991    12165120 :          speccomb(lay) = colh2o(lay) + rat_h2och4(lay)*colch4(lay)
    2992    12165120 :          specparm(lay) = colh2o(lay)/speccomb(lay)
    2993    12165120 :          if (specparm(lay) .ge. oneminus) specparm(lay) = oneminus
    2994    12165120 :          specmult(lay) = 8._r8*(specparm(lay))
    2995    12165120 :          js(lay) = 1 + int(specmult(lay))
    2996    12165120 :          fs(lay) = mod(specmult(lay),1.0_r8)
    2997             : 
    2998    12165120 :          speccomb1(lay) = colh2o(lay) + rat_h2och4_1(lay)*colch4(lay)
    2999    12165120 :          specparm1(lay) = colh2o(lay)/speccomb1(lay)
    3000    12165120 :          if (specparm1(lay) .ge. oneminus) specparm1(lay) = oneminus
    3001    12165120 :          specmult1(lay) = 8._r8*(specparm1(lay))
    3002    12165120 :          js1(lay) = 1 + int(specmult1(lay))
    3003    12165120 :          fs1(lay) = mod(specmult1(lay),1.0_r8)
    3004             : 
    3005    12165120 :          speccomb_planck(lay) = colh2o(lay)+refrat_planck_a*colch4(lay)
    3006    12165120 :          specparm_planck(lay) = colh2o(lay)/speccomb_planck(lay)
    3007    12165120 :          if (specparm_planck(lay) .ge. oneminus) specparm_planck(lay)=oneminus
    3008    12165120 :          specmult_planck(lay) = 8._r8*specparm_planck(lay)
    3009    12165120 :          jpl(lay)= 1 + int(specmult_planck(lay))
    3010    12165120 :          fpl(lay) = mod(specmult_planck(lay),1.0_r8)
    3011             : 
    3012    12165120 :          ind0(lay) = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(16) + js(lay)
    3013    12718080 :          ind1(lay) = (jp(lay)*5+(jt1(lay)-1))*nspa(16) + js1(lay)
    3014             :       enddo
    3015    12718080 :       do lay = 1, laytrop
    3016    12718080 :          if (specparm(lay) .lt. 0.125_r8) then
    3017     2156503 :             p = fs(lay) - 1
    3018     2156503 :             p4 = p**4
    3019     2156503 :             fk0 = p4
    3020     2156503 :             fk1 = 1 - p - 2.0_r8*p4
    3021     2156503 :             fk2 = p + p4
    3022     2156503 :             fac000(lay) = fk0*fac00(lay)
    3023     2156503 :             fac100(lay) = fk1*fac00(lay)
    3024     2156503 :             fac200(lay) = fk2*fac00(lay)
    3025     2156503 :             fac010(lay) = fk0*fac10(lay)
    3026     2156503 :             fac110(lay) = fk1*fac10(lay)
    3027     2156503 :             fac210(lay) = fk2*fac10(lay)
    3028    10008617 :          else if (specparm(lay) .gt. 0.875_r8) then
    3029       29999 :             p = -fs(lay) 
    3030       29999 :             p4 = p**4
    3031       29999 :             fk0 = p4
    3032       29999 :             fk1 = 1 - p - 2.0_r8*p4
    3033       29999 :             fk2 = p + p4
    3034       29999 :             fac000(lay) = fk0*fac00(lay)
    3035       29999 :             fac100(lay) = fk1*fac00(lay)
    3036       29999 :             fac200(lay) = fk2*fac00(lay)
    3037       29999 :             fac010(lay) = fk0*fac10(lay)
    3038       29999 :             fac110(lay) = fk1*fac10(lay)
    3039       29999 :             fac210(lay) = fk2*fac10(lay)
    3040             :          else
    3041     9978618 :             fac000(lay) = (1._r8 - fs(lay)) * fac00(lay)
    3042     9978618 :             fac010(lay) = (1._r8 - fs(lay)) * fac10(lay)
    3043     9978618 :             fac100(lay) = fs(lay) * fac00(lay)
    3044     9978618 :             fac110(lay) = fs(lay) * fac10(lay)
    3045             :          endif
    3046             :       enddo
    3047    12718080 :       do lay = 1, laytrop
    3048             : 
    3049    12718080 :          if (specparm1(lay) .lt. 0.125_r8) then
    3050      845664 :             p = fs1(lay) - 1
    3051      845664 :             p4 = p**4
    3052      845664 :             fk0 = p4
    3053      845664 :             fk1 = 1 - p - 2.0_r8*p4
    3054      845664 :             fk2 = p + p4
    3055      845664 :             fac001(lay) = fk0*fac01(lay)
    3056      845664 :             fac101(lay) = fk1*fac01(lay)
    3057      845664 :             fac201(lay) = fk2*fac01(lay)
    3058      845664 :             fac011(lay) = fk0*fac11(lay)
    3059      845664 :             fac111(lay) = fk1*fac11(lay)
    3060      845664 :             fac211(lay) = fk2*fac11(lay)
    3061    11319456 :          else if (specparm1(lay) .gt. 0.875_r8) then
    3062      365328 :             p = -fs1(lay) 
    3063      365328 :             p4 = p**4
    3064      365328 :             fk0 = p4
    3065      365328 :             fk1 = 1 - p - 2.0_r8*p4
    3066      365328 :             fk2 = p + p4
    3067      365328 :             fac001(lay) = fk0*fac01(lay)
    3068      365328 :             fac101(lay) = fk1*fac01(lay)
    3069      365328 :             fac201(lay) = fk2*fac01(lay)
    3070      365328 :             fac011(lay) = fk0*fac11(lay)
    3071      365328 :             fac111(lay) = fk1*fac11(lay)
    3072      365328 :             fac211(lay) = fk2*fac11(lay)
    3073             :          else
    3074    10954128 :             fac001(lay) = (1._r8 - fs1(lay)) * fac01(lay)
    3075    10954128 :             fac011(lay) = (1._r8 - fs1(lay)) * fac11(lay)
    3076    10954128 :             fac101(lay) = fs1(lay) * fac01(lay)
    3077    10954128 :             fac111(lay) = fs1(lay) * fac11(lay)
    3078             :          endif
    3079             :       enddo
    3080    12718080 :       do lay = 1, laytrop
    3081             : 
    3082    37048320 :          do ig = 1, ng16
    3083    72990720 :             tauself = selffac(lay)* (selfref(indself(lay),ig) + selffrac(lay) * &
    3084    97320960 :                  (selfref(indself(lay)+1,ig) - selfref(indself(lay),ig)))
    3085    48660480 :             taufor =  forfac(lay) * (forref(indfor(lay),ig) + forfrac(lay) * &
    3086    48660480 :                  (forref(indfor(lay)+1,ig) - forref(indfor(lay),ig))) 
    3087             : 
    3088    24330240 :             if (specparm(lay) .lt. 0.125_r8) then
    3089             :                tau_major = speccomb(lay) * &
    3090     4313006 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3091     4313006 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3092     4313006 :                     fac200(lay) * absa(ind0(lay)+2,ig) + &
    3093     4313006 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3094     4313006 :                     fac110(lay) * absa(ind0(lay)+10,ig) + &
    3095    25878036 :                     fac210(lay) * absa(ind0(lay)+11,ig))
    3096    20017234 :             else if (specparm(lay) .gt. 0.875_r8) then
    3097             :                tau_major = speccomb(lay) * &
    3098       59998 :                     (fac200(lay) * absa(ind0(lay)-1,ig) + &
    3099       59998 :                     fac100(lay) * absa(ind0(lay),ig) + &
    3100       59998 :                     fac000(lay) * absa(ind0(lay)+1,ig) + &
    3101       59998 :                     fac210(lay) * absa(ind0(lay)+8,ig) + &
    3102       59998 :                     fac110(lay) * absa(ind0(lay)+9,ig) + &
    3103      359988 :                     fac010(lay) * absa(ind0(lay)+10,ig))
    3104             :             else
    3105             :                tau_major = speccomb(lay) * &
    3106    19957236 :                     (fac000(lay) * absa(ind0(lay),ig) + &
    3107    19957236 :                     fac100(lay) * absa(ind0(lay)+1,ig) + &
    3108    19957236 :                     fac010(lay) * absa(ind0(lay)+9,ig) + &
    3109    79828944 :                     fac110(lay) * absa(ind0(lay)+10,ig))
    3110             :             endif
    3111             : 
    3112    24330240 :             if (specparm1(lay) .lt. 0.125_r8) then
    3113             :                tau_major1 = speccomb1(lay) * &
    3114     1691328 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3115     1691328 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3116     1691328 :                     fac201(lay) * absa(ind1(lay)+2,ig) + &
    3117     1691328 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3118     1691328 :                     fac111(lay) * absa(ind1(lay)+10,ig) + &
    3119    10147968 :                     fac211(lay) * absa(ind1(lay)+11,ig))
    3120    22638912 :             else if (specparm1(lay) .gt. 0.875_r8) then
    3121             :                tau_major1 = speccomb1(lay) * &
    3122      730656 :                     (fac201(lay) * absa(ind1(lay)-1,ig) + &
    3123      730656 :                     fac101(lay) * absa(ind1(lay),ig) + &
    3124      730656 :                     fac001(lay) * absa(ind1(lay)+1,ig) + &
    3125      730656 :                     fac211(lay) * absa(ind1(lay)+8,ig) + &
    3126      730656 :                     fac111(lay) * absa(ind1(lay)+9,ig) + &
    3127     4383936 :                     fac011(lay) * absa(ind1(lay)+10,ig))
    3128             :             else
    3129             :                tau_major1 = speccomb1(lay) * &
    3130    21908256 :                     (fac001(lay) * absa(ind1(lay),ig) + &
    3131    21908256 :                     fac101(lay) * absa(ind1(lay)+1,ig) + &
    3132    21908256 :                     fac011(lay) * absa(ind1(lay)+9,ig) + &
    3133    87633024 :                     fac111(lay) * absa(ind1(lay)+10,ig))
    3134             :             endif
    3135             : 
    3136    24330240 :             taug(lay,ngs15+ig) = tau_major + tau_major1 &
    3137    24330240 :                  + tauself + taufor
    3138    48660480 :             fracs(lay,ngs15+ig) = fracrefa(ig,jpl(lay)) + fpl(lay) * &
    3139    60825600 :                  (fracrefa(ig,jpl(lay)+1)-fracrefa(ig,jpl(lay)))
    3140             :          enddo
    3141             :       enddo
    3142             : 
    3143             : ! Upper atmosphere loop
    3144    18800640 :       do lay = laytrop+1, nlayers
    3145    18247680 :          ind0(lay) = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(16) + 1
    3146    18800640 :          ind1(lay) = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(16) + 1
    3147             :       enddo
    3148    18800640 :       do lay = laytrop+1, nlayers
    3149    55296000 :          do ig = 1, ng16
    3150    72990720 :             taug(lay,ngs15+ig) = colch4(lay) * &
    3151   109486080 :                  (fac00(lay) * absb(ind0(lay),ig) + &
    3152    72990720 :                  fac10(lay) * absb(ind0(lay)+1,ig) + &
    3153    72990720 :                  fac01(lay) * absb(ind1(lay),ig) + &
    3154   255467520 :                  fac11(lay) * absb(ind1(lay)+1,ig))
    3155    54743040 :             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