!     path:      $Source$
!     author:    $Author: mike $
!     revision:  $Revision: 11661 $
!     created:   $Date: 2009-05-22 18:22:22 -0400 (Fri, 22 May 2009) $

      module rrtmg_sw_taumol
!  --------------------------------------------------------------------------
! |                                                                          |
! |  Copyright 2002-2009, Atmospheric & Environmental Research, Inc. (AER).  |
! |  This software may be used, copied, or redistributed as long as it is    |
! |  not sold and this copyright notice is reproduced on each copy made.     |
! |  This model is provided as is without any express or implied warranties. |
! |                       (http://www.rtweb.aer.com/)                        |
! |                                                                          |
!  --------------------------------------------------------------------------

! ------- Modules -------

      use parkind,  only : im => kind_im, rb => kind_rb
      use parrrsw,  only : ng24, ng25
      use rrsw_con, only:  oneminus
      use rrsw_wvn, only: nspa, nspb
      use rrsw_vsn, only: hvrtau, hnamtau
      use rrsw_fastj_cmn
      use fjx_cmn_mod, only: W_r

      implicit none

      contains

!----------------------------------------------------------------------------
      subroutine taumol_sw(nlayers, &
                           colh2o, colco2, colch4, colo2, colo3, colmol, &
                           laytrop, jp, jt, jt1, &
                           fac00, fac01, fac10, fac11, &
                           selffac, selffrac, indself, forfac, forfrac, indfor, &
                           taug_rrtmg)
!----------------------------------------------------------------------------
! ******************************************************************************
! *                                                                            *
! *                 Optical depths developed for the                           *
! *                                                                            *
! *               RAPID RADIATIVE TRANSFER MODEL (RRTM)                        *
! *                                                                            *
! *                                                                            *
! *           ATMOSPHERIC AND ENVIRONMENTAL RESEARCH, INC.                     *
! *                       131 HARTWELL AVENUE                                  *
! *                       LEXINGTON, MA 02421                                  *
! *                                                                            *
! *                                                                            *
! *                          ELI J. MLAWER                                     *
! *                        JENNIFER DELAMERE                                   *
! *                        STEVEN J. TAUBMAN                                   *
! *                        SHEPARD A. CLOUGH                                   *
! *                                                                            *
! *                                                                            *
! *                                                                            *
! *                                                                            *
! *                      email:  mlawer@aer.com                                *
! *                      email:  jdelamer@aer.com                              *
! *                                                                            *
! *       The authors wish to acknowledge the contributions of the             *
! *       following people:  Patrick D. Brown, Michael J. Iacono,              *
! *       Ronald E. Farren, Luke Chen, Robert Bergstrom.                       *
! *                                                                            *
! ******************************************************************************
! *    TAUMOL                                                                  *
! *                                                                            *
! *    This file contains the subroutines TAUGBn (where n goes from            *
! *    1 to 28).  TAUGBn calculates the optical depths and Planck fractions    *
! *    per g-value and layer for band n.                                       *
! *                                                                            *
! * Output:  optical depths (unitless)                                         *
! *          fractions needed to compute Planck functions at every layer       *
! *              and g-value                                                   *
! *                                                                            *
! *    COMMON /TAUGCOM/  TAUG(MXLAY,MG)                                        *
! *    COMMON /PLANKG/   FRACS(MXLAY,MG)                                       *
! *                                                                            *
! * Input                                                                      *
! *                                                                            *
! *    PARAMETER (MG=16, MXLAY=203, NBANDS=14)                                 *
! *                                                                            *
! *    COMMON /FEATURES/ NG(NBANDS),NSPA(NBANDS),NSPB(NBANDS)                  *
! *    COMMON /PRECISE/  ONEMINUS                                              *
! *    COMMON /PROFILE/  NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),                    *
! *   &                  PZ(0:MXLAY),TZ(0:MXLAY),TBOUND                        *
! *    COMMON /PROFDATA/ LAYTROP,LAYSWTCH,LAYLOW,                              *
! *   &                  COLH2O(MXLAY),COLCO2(MXLAY),                          *
! *   &                  COLO3(MXLAY),COLN2O(MXLAY),COLCH4(MXLAY),             *
! *   &                  COLO2(MXLAY),CO2MULT(MXLAY)                           *
! *    COMMON /INTFAC/   FAC00(MXLAY),FAC01(MXLAY),                            *
! *   &                  FAC10(MXLAY),FAC11(MXLAY)                             *
! *    COMMON /INTIND/   JP(MXLAY),JT(MXLAY),JT1(MXLAY)                        *
! *    COMMON /SELF/     SELFFAC(MXLAY), SELFFRAC(MXLAY), INDSELF(MXLAY)       *
! *                                                                            *
! *    Description:                                                            *
! *    NG(IBAND) - number of g-values in band IBAND                            *
! *    NSPA(IBAND) - for the lower atmosphere, the number of reference         *
! *                  atmospheres that are stored for band IBAND per            *
! *                  pressure level and temperature.  Each of these            *
! *                  atmospheres has different relative amounts of the         *
! *                  key species for the band (i.e. different binary           *
! *                  species parameters).                                      *
! *    NSPB(IBAND) - same for upper atmosphere                                 *
! *    ONEMINUS - since problems are caused in some cases by interpolation     *
! *               parameters equal to or greater than 1, for these cases       *
! *               these parameters are set to this value, slightly < 1.        *
! *    PAVEL - layer pressures (mb)                                            *
! *    TAVEL - layer temperatures (degrees K)                                  *
! *    PZ - level pressures (mb)                                               *
! *    TZ - level temperatures (degrees K)                                     *
! *    LAYTROP - layer at which switch is made from one combination of         *
! *              key species to another                                        *
! *    COLH2O, COLCO2, COLO3, COLN2O, COLCH4 - column amounts of water         *
! *              vapor,carbon dioxide, ozone, nitrous ozide, methane,          *
! *              respectively (molecules/cm**2)                                *
! *    CO2MULT - for bands in which carbon dioxide is implemented as a         *
! *              trace species, this is the factor used to multiply the        *
! *              band's average CO2 absorption coefficient to get the added    *
! *              contribution to the optical depth relative to 355 ppm.        *
! *    FACij(LAY) - for layer LAY, these are factors that are needed to        *
! *                 compute the interpolation factors that multiply the        *
! *                 appropriate reference k-values.  A value of 0 (1) for      *
! *                 i,j indicates that the corresponding factor multiplies     *
! *                 reference k-value for the lower (higher) of the two        *
! *                 appropriate temperatures, and altitudes, respectively.     *
! *    JP - the index of the lower (in altitude) of the two appropriate        *
! *         reference pressure levels needed for interpolation                 *
! *    JT, JT1 - the indices of the lower of the two appropriate reference     *
! *              temperatures needed for interpolation (for pressure           *
! *              levels JP and JP+1, respectively)                             *
! *    SELFFAC - scale factor needed to water vapor self-continuum, equals     *
! *              (water vapor density)/(atmospheric density at 296K and        *
! *              1013 mb)                                                      *
! *    SELFFRAC - factor needed for temperature interpolation of reference     *
! *               water vapor self-continuum data                              *
! *    INDSELF - index of the lower of the two appropriate reference           *
! *              temperatures needed for the self-continuum interpolation      *
! *                                                                            *
! * Data input                                                                 *
! *    COMMON /Kn/ KA(NSPA(n),5,13,MG), KB(NSPB(n),5,13:59,MG), SELFREF(10,MG) *
! *       (note:  n is the band number)                                        *
! *                                                                            *
! *    Description:                                                            *
! *    KA - k-values for low reference atmospheres (no water vapor             *
! *         self-continuum) (units: cm**2/molecule)                            *
! *    KB - k-values for high reference atmospheres (all sources)              *
! *         (units: cm**2/molecule)                                            *
! *    SELFREF - k-values for water vapor self-continuum for reference         *
! *              atmospheres (used below LAYTROP)                              *
! *              (units: cm**2/molecule)                                       *
! *                                                                            *
! *    DIMENSION ABSA(65*NSPA(n),MG), ABSB(235*NSPB(n),MG)                     *
! *    EQUIVALENCE (KA,ABSA),(KB,ABSB)                                         *
! *                                                                            *
! *****************************************************************************
!
! Modifications
!
! Revised: Adapted to F90 coding, J.-J.Morcrette, ECMWF, Feb 2003
! Revised: Modified for g-point reduction, MJIacono, AER, Dec 2003
! Revised: Reformatted for consistency with rrtmg_lw, MJIacono, AER, Jul 2006
!
! ------- Declarations -------

! ----- Input -----
      integer(kind=im), intent(in) :: nlayers            ! total number of layers

      integer(kind=im), intent(in) :: laytrop            ! tropopause layer index
      integer(kind=im), intent(in) :: jp(:)              ! 
                                                         !   Dimensions: (nlayers)
      integer(kind=im), intent(in) :: jt(:)              !
                                                         !   Dimensions: (nlayers)
      integer(kind=im), intent(in) :: jt1(:)             !
                                                         !   Dimensions: (nlayers)

      real(kind=rb), intent(in) :: colh2o(:)             ! column amount (h2o)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: colco2(:)             ! column amount (co2)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: colo3(:)              ! column amount (o3)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: colch4(:)             ! column amount (ch4)
                                                         !   Dimensions: (nlayers)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: colo2(:)              ! column amount (o2)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: colmol(:)             ! 
                                                         !   Dimensions: (nlayers)

      integer(kind=im), intent(in) :: indself(:)    
                                                         !   Dimensions: (nlayers)
      integer(kind=im), intent(in) :: indfor(:)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: selffac(:)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: selffrac(:)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: forfac(:)
                                                         !   Dimensions: (nlayers)
      real(kind=rb), intent(in) :: forfrac(:)
                                                         !   Dimensions: (nlayers)

      real(kind=rb), intent(in) :: &                     !
                       fac00(:), fac01(:), &             !   Dimensions: (nlayers)
                       fac10(:), fac11(:) 

! ----- Output -----

                                                         !   Dimensions: (ngptsw)
      real(kind=rb), intent(out) :: taug_rrtmg(:, 0:)            ! gaseous optical depth 

!      real(kind=rb), intent(out) :: sfluxzen(:)          ! solar source function
                                                         !   Dimensions: (nlayers,ngptsw)
!      real(kind=rb), intent(out) :: taur_rrtmg(:,:)            ! Rayleigh 
                                                         !   Dimensions: (nlayers,ngptsw)
!      real(kind=rb), intent(out) :: ssa(:,:)            ! single scattering albedo (inactive)
                                                         !   Dimensions: (nlayers,ngptsw)
      
!local 
      real(kind=rb) :: taug24(nlayers, ng24), taug25(nlayers, ng25)
!      real(kind=rb) :: taur24(nlayers, ng24), taur25(nlayers, ng25)
      real(kind=rb) :: sflx24(ng24),  sflx25(ng25)

      integer(kind=im) :: kout 


      hvrtau = '$Revision: 11661 $'

      oneminus= 1.0_rb - 1.e-6_rb
! Calculate gaseous optical depth and planck fractions for each spectral band.
      taug_rrtmg(:,:)=0.d0
!
! reorder from short wave to long wave length
!
      call taumol24(taug24, sflx24)
      call taumol25(taug25, sflx25)
      call taumol_c(taug24,sflx24, taug25,sflx25, kout)
      call taumol23(kout)
      call taumol22(kout)
      call taumol21(kout)
      call taumol20(kout)
      call taumol19(kout)
      call taumol18(kout)
      call taumol17(kout)
      call taumol16(kout)
      call taumol29(kout)
      if (kout .ne. W_r) print*,'kout=', kout, 'W_r=', W_r
!-------------
      contains
!-------------


      subroutine taumol_c(taug24, sflux24, taug25, sflux25, kout) 
        
      use parrrsw, only : ng24, ng25
      
      real(kind=rb), intent(in) :: taug24(nlayers, ng24), sflux24(ng24)
      real(kind=rb), intent(in) :: taug25(nlayers, ng25), sflux25(ng25)    
!      
      integer(kind=im), intent(out) :: kout 
!
      real(kind=rb):: wtaug, sf
      integer(kind=im)::  I, L
! disable this one
!      DO L=1,nlayers
!         wtaug=0.d0
!         wtaur=0.d0
!         sf= 0.d0
!         DO I=1,5
!            wtaug =  taug25(L,I)* sflux25(I)+ wtaug
!             sf   =  sflux25(I)+ sf
!         ENDDO
!         DO I=1,4
!            wtaug =  taug24(L,I)* sflux24(I)+ wtaug
!             sf   =  sflux24(I)+ sf
!         ENDDO
!         taug_rrtmg(L,1)= wtaug/sf
!         taur_rrtmg(L,1)= wtaur/sf
!      ENDDO


      DO L=1, nlayers
         wtaug=0.d0
         sf=0.d0
         wtaug = taug25(L,5)* sflux25(5) + &
                 taug25(L,4)* sflux25(4) + &
                 taug25(L,3)* sflux25(3) + &
                 taug25(L,2)* sflux25(2) + &
                 taug25(L,1)* sflux25(1) + &
                 taug24(L,4)* sflux24(4) + &
                 taug24(L,3)* sflux24(3) + &
                 taug24(L,2)* sflux24(2) + &
                 taug24(L,1)* sflux24(1)
!   
         sf=                  sflux25(5) + &      
                              sflux25(4) + &
                              sflux25(3) + &
                              sflux25(2) + &
                              sflux25(1) + &
                              sflux24(4) + &
                              sflux24(3) + &
                              sflux24(2) + &
                              sflux24(1)   

         taug_rrtmg(L,0) = wtaug/sf* 0.9d0! this will add onto SolarJ bin-18 (under O3 Chappius band)
      ENDDO
      DO L=1,nlayers
         wtaug=0.d0
         sf= 0.d0
         wtaug  =   taug25(L,6)* sflux25(6)+ &
                    taug24(L,5)* sflux24(5)
         sf   =     sflux25(6)+ sflux24(5)
         taug_rrtmg(L,1)= wtaug/sf
      ENDDO
      DO L=1,nlayers
         taug_rrtmg(L,2)= taug24(L,6)
         taug_rrtmg(L,3)= taug24(L,7)
         taug_rrtmg(L,4)= taug24(L,8)
      ENDDO
      kout= 4       
      end subroutine taumol_c          


!----------------------------------------------------------------------------
      subroutine taumol16(kout)
!----------------------------------------------------------------------------
!
!     band 16:  2600-3250 cm-1 (low - h2o,ch4; high - ch4)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng16, ngs17
      use rrsw_kg16, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
      integer(kind=im), intent(inout):: kout

! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray, strrat1

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      strrat1 = 252.131_rb
      layreffr = 18

! Lower atmosphere loop
      do lay = 1, laytrop
         speccomb = colh2o(lay) + strrat1*colch4(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(16) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(16) + js
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl

         do ig = 1, ng16
            taug_rrtmg(lay,ngs17+ig) = speccomb * &
                (fac000 * absa(ind0   ,ig) + &
                 fac100 * absa(ind0 +1,ig) + &
                 fac010 * absa(ind0 +9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1   ,ig) + &
                 fac101 * absa(ind1 +1,ig) + &
                 fac011 * absa(ind1 +9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
                 colh2o(lay) * &
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) 
!            ssa(lay,ngs17+ ig) = tauray/taug_rrtmg(lay,ig)
!            taur_rrtmg(lay,ngs17+ig) = tauray
         enddo
      enddo

      laysolfr = nlayers

! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
            laysolfr = lay
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(16) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(16) + 1
         tauray = colmol(lay) * rayl
         do ig = 1, ng16
            taug_rrtmg(lay,ngs17+ig) = colch4(lay) * &
                (fac00(lay) * absb(ind0  ,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1  ,ig) + &
                 fac11(lay) * absb(ind1+1,ig)) 
!             ssa(lay,ngs17+ig) = tauray/taug_rrtmg(lay,ig)
!             if (lay .eq. laysolfr) sfluxzen(ngs17+ig) = sfluxref(ig) 
!             taur_rrtmg(lay,ngs17+ ig) = tauray  
         enddo
      enddo
      kout= kout + ng16
      end subroutine taumol16

!----------------------------------------------------------------------------
      subroutine taumol17(kout)
!----------------------------------------------------------------------------
!
!     band 17:  3250-4000 cm-1 (low - h2o,co2; high - h2o,co2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng17, ngs18
      use rrsw_kg17, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
      integer(kind=im), intent(inout):: kout
! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray, strrat

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      strrat = 0.364641_rb
      layreffr = 30

! Lower atmosphere loop
      do lay = 1, laytrop
         speccomb = colh2o(lay) + strrat*colco2(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(17) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(17) + js
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl
         do ig = 1, ng17
            taug_rrtmg(lay,ngs18+ig) = speccomb * &
                 (fac000 * absa(ind0,ig) + &
                 fac100 * absa(ind0+1,ig) + &
                 fac010 * absa(ind0+9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1,ig) + &
                 fac101 * absa(ind1+1,ig) + &
                 fac011 * absa(ind1+9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
                 colh2o(lay) * &
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) 
!            ssa(lay,ngs18+ig) = tauray/taug_rrtmg(lay,ngs18+ig)
!            taur_rrtmg(lay,ngs18+ig) = tauray
         enddo
      enddo
      laysolfr = nlayers
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
            laysolfr = lay
         speccomb = colh2o(lay) + strrat*colco2(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 4._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1)) *nspb(17) + js
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(17) + js
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl
         do ig = 1, ng17
            taug_rrtmg(lay,ngs18+ig) = speccomb * &
                (fac000 * absb(ind0,  ig) + &
                 fac100 * absb(ind0+1,ig) + &
                 fac010 * absb(ind0+5,ig) + &
                 fac110 * absb(ind0+6,ig) + &
                 fac001 * absb(ind1,ig) + &
                 fac101 * absb(ind1+1,ig) + &
                 fac011 * absb(ind1+5,ig) + &
                 fac111 * absb(ind1+6,ig)) + &
                 colh2o(lay) * &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig))) 
!            ssa(lay,ngs18+ig) = tauray/taug_rrtmg(lay,ngs18+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs18+ig) = sfluxref(ig,js) &
!                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taur_rrtmg(lay,ngs18+ig) = tauray
         enddo
      enddo
      kout = kout + ng17
      end subroutine taumol17

!----------------------------------------------------------------------------
      subroutine taumol18(kout)
!----------------------------------------------------------------------------
!
!     band 18:  4000-4650 cm-1 (low - h2o,ch4; high - ch4)
!
!----------------------------------------------------------------------------
! ------- Modules -------

      use parrrsw, only : ng18, ngs19
      use rrsw_kg18, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
      integer(kind= im), intent(inout):: kout      
! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray, strrat

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      strrat = 38.9589_rb
      layreffr = 6
      laysolfr = laytrop
      
! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
            laysolfr = min(lay+1,laytrop)
         speccomb = colh2o(lay) + strrat*colch4(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(18) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(18) + js
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl

         do ig = 1, ng18
            taug_rrtmg(lay,ngs19+ig) = speccomb * &
                (fac000 * absa(ind0,ig) + &
                 fac100 * absa(ind0+1,ig) + &
                 fac010 * absa(ind0+9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1,ig) + &
                 fac101 * absa(ind1+1,ig) + &
                 fac011 * absa(ind1+9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
                 colh2o(lay) * &
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) 
!            ssa(lay,ngs19+ig) = tauray/taug_rrtmg(lay,ngs19+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs19+ig) = sfluxref(ig,js) &
!                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taur_rrtmg(lay,ngs19+ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(18) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(18) + 1
!         tauray = colmol(lay) * rayl
         do ig = 1, ng18
            taug_rrtmg(lay,ngs19+ig) = colch4(lay) * &
                (fac00(lay) * absb(ind0,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1,ig) + &	  
                 fac11(lay) * absb(ind1+1,ig)) 
 !           ssa(lay,ngs19+ig) = tauray/taug_rrtmg(lay,ngs19+ig)
 !           taur_rrtmg(lay,ngs19+ig) = tauray
         enddo
       enddo
       kout= kout + ng18
       end subroutine taumol18

!----------------------------------------------------------------------------
      subroutine taumol19(kout)
!----------------------------------------------------------------------------
!
!     band 19:  4650-5150 cm-1 (low - h2o,co2; high - co2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

        use parrrsw, only : ng19, ngs20
        use rrsw_kg19, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
        integer(kind=im), intent(inout):: kout
! Local

        integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
             layreffr
        real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
             fac110, fac111, fs, speccomb, specmult, specparm, &
             tauray, strrat
! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  
        strrat = 5.49281_rb
        layreffr = 3
        laysolfr = laytrop
! Lower atmosphere loop      
        do lay = 1, laytrop
           if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
                laysolfr = min(lay+1,laytrop)
           speccomb = colh2o(lay) + strrat*colco2(lay)
           specparm = colh2o(lay)/speccomb 
           if (specparm .ge. oneminus) specparm = oneminus
           specmult = 8._rb*(specparm)
           js = 1 + int(specmult)
           fs = mod(specmult, 1._rb )
           fac000 = (1._rb - fs) * fac00(lay)
           fac010 = (1._rb - fs) * fac10(lay)
           fac100 = fs * fac00(lay)
           fac110 = fs * fac10(lay)
           fac001 = (1._rb - fs) * fac01(lay)
           fac011 = (1._rb - fs) * fac11(lay)
           fac101 = fs * fac01(lay)
           fac111 = fs * fac11(lay)
           ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(19) + js
           ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(19) + js
           inds = indself(lay)
           indf = indfor(lay)
!           tauray = colmol(lay) * rayl
           do ig = 1 , ng19
              taug_rrtmg(lay,ngs20+ig) = speccomb * &
                   (fac000 * absa(ind0,ig) + &
                   fac100 * absa(ind0+1,ig) + &
                   fac010 * absa(ind0+9,ig) + &
                   fac110 * absa(ind0+10,ig) + &
                   fac001 * absa(ind1,ig) + &
                   fac101 * absa(ind1+1,ig) + &
                   fac011 * absa(ind1+9,ig) + &
                   fac111 * absa(ind1+10,ig)) + &
                   colh2o(lay) * &
                   (selffac(lay) * (selfref(inds,ig) + &
                   selffrac(lay) * &
                   (selfref(inds+1,ig) - selfref(inds,ig))) + & 
                   forfac(lay) * (forref(indf,ig) + &
                   forfrac(lay) * &
                   (forref(indf+1,ig) - forref(indf,ig)))) 
!            ssa(lay,ngs20+ig) = tauray/taug_rrtmg(lay,ngs20+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs20 +ig) = sfluxref(ig,js) &
!                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taur_rrtmg(lay,ngs20+ig) = tauray   
           enddo
        enddo
! Upper atmosphere loop
        do lay = laytrop+1, nlayers
           ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(19) + 1
           ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(19) + 1
!           tauray = colmol(lay) * rayl
           do ig = 1 , ng19
              taug_rrtmg(lay,ngs20+ig) = colco2(lay) * &
                   (fac00(lay) * absb(ind0,ig) + &
                   fac10(lay) * absb(ind0+1,ig) + &
                   fac01(lay) * absb(ind1,ig) + &
                   fac11(lay) * absb(ind1+1,ig)) 
!            ssa(lay,ngs20+ig) = tauray/taug_rrtmg(lay,ngs20+ig) 
!            taur_rrtmg(lay,ngs20+ig) = tauray   
           enddo
        enddo
        kout= kout+ ng19
      end subroutine taumol19
!----------------------------------------------------------------------------
      subroutine taumol20(kout)
!----------------------------------------------------------------------------
!
!     band 20:  5150-6150 cm-1 (low - h2o; high - h2o)
!
!----------------------------------------------------------------------------

! ------- Modules -------

        use parrrsw,   only : ng20, ngs21
        use rrsw_kg20, only : absa, ka, absb, kb, forref, selfref, &
             sfluxref, absch4, rayl

        implicit none

! ------- Declarations -------
        integer(kind=im), intent(inout):: kout 
! Local

        integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
             layreffr
        real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
             fac110, fac111, fs, speccomb, specmult, specparm, &
             tauray

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

        layreffr = 3
        laysolfr = laytrop

! Lower atmosphere loop
        do lay = 1, laytrop
           if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
                laysolfr = min(lay+1,laytrop)
           ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(20) + 1
           ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(20) + 1
           inds = indself(lay)
           indf = indfor(lay)
!           tauray = colmol(lay) * rayl
           do ig = 1, ng20
              taug_rrtmg(lay,ngs21+ig) = colh2o(lay) * &
                   ((fac00(lay) * absa(ind0,ig) + &
                   fac10(lay) * absa(ind0+1,ig) + &
                   fac01(lay) * absa(ind1,ig) + &
                   fac11(lay) * absa(ind1+1,ig)) + &
                   selffac(lay) * (selfref(inds,ig) + & 
                   selffrac(lay) * &
                   (selfref(inds+1,ig) - selfref(inds,ig))) + &
                   forfac(lay) * (forref(indf,ig) + &
                   forfrac(lay) * &
                   (forref(indf+1,ig) - forref(indf,ig)))) &
                   + colch4(lay) * absch4(ig)
!           ssa(lay,ngs21+ig) = tauray/taug_rrtmg(lay,ngs21+ig)
!            taur_rrtmg(lay,ngs21+ig) = tauray 
!            if (lay .eq. laysolfr) sfluxzen(ngs21+ig) = sfluxref(ig) 
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(20) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(20) + 1
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl
         do ig = 1, ng20
            taug_rrtmg(lay,ngs21+ig) = colh2o(lay) * &
                 (fac00(lay) * absb(ind0,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1,ig) + &
                 fac11(lay) * absb(ind1+1,ig) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) + &
                 colch4(lay) * absch4(ig)
!            ssa(lay,ngs21+ig) = tauray/taug_rrtmg(lay,ngs21+ig)
!            taur_rrtmg(lay,ngs21+ig) = tauray 
         enddo
      enddo
      kout = kout + ng20
    end subroutine taumol20

!----------------------------------------------------------------------------
    subroutine taumol21(kout)
!----------------------------------------------------------------------------
!
!     band 21:  6150-7700 cm-1 (low - h2o,co2; high - h2o,co2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng21, ngs22
      use rrsw_kg21, only : absa, ka, absb, kb, forref, selfref, &
           sfluxref, rayl

! ------- Declarations -------
      integer(kind=im), intent(inout):: kout

! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
           layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
           fac110, fac111, fs, speccomb, specmult, specparm, &
           tauray, strrat
! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      strrat = 0.0045321_rb
      layreffr = 8
      laysolfr = laytrop      
! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
              laysolfr = min(lay+1,laytrop)
         speccomb = colh2o(lay) + strrat*colco2(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(21) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(21) + js
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl

         do ig = 1, ng21
            taug_rrtmg(lay,ngs22+ig) = speccomb * &
                (fac000 * absa(ind0,ig) + &
                 fac100 * absa(ind0+1,ig) + &
                 fac010 * absa(ind0+9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1,ig) + &
                 fac101 * absa(ind1+1,ig) + &
                 fac011 * absa(ind1+9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
                 colh2o(lay) * &
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig))))
!            ssa(lay,ngs22+ig) = tauray/taug_rrtmg(lay,ngs22+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs22+ig) = sfluxref(ig,js) &
!                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taur_rrtmg(lay,ngs22+ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         speccomb = colh2o(lay) + strrat*colco2(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 4._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(21) + js
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(21) + js
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl
         do ig = 1, ng21
            taug_rrtmg(lay,ngs22+ig) = speccomb * &
                (fac000 * absb(ind0,ig) + &
                 fac100 * absb(ind0+1,ig) + &
                 fac010 * absb(ind0+5,ig) + &
                 fac110 * absb(ind0+6,ig) + &
                 fac001 * absb(ind1,ig) + &
                 fac101 * absb(ind1+1,ig) + &
                 fac011 * absb(ind1+5,ig) + &
                 fac111 * absb(ind1+6,ig)) + &
                 colh2o(lay) * &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))
!            ssa(lay,ngs22+ig) = tauray/taug_rrtmg(lay,ngs22+ig)
!            taur_rrtmg(lay,ngs22+ig) = tauray
         enddo
      enddo
      kout=kout + ng21
      end subroutine taumol21

!----------------------------------------------------------------------------
      subroutine taumol22(kout)
!----------------------------------------------------------------------------
!
!     band 22:  7700-8050 cm-1 (low - h2o,o2; high - o2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng22, ngs23
      use rrsw_kg22, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
      integer(kind=im),intent(inout):: kout

! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray, o2adj, o2cont, strrat

! The following factor is the ratio of total O2 band intensity (lines 
! and Mate continuum) to O2 band intensity (line only).  It is needed
! to adjust the optical depths since the k's include only lines.
      o2adj = 1.6_rb
      
! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      strrat = 0.022708_rb
      layreffr = 2
      laysolfr = laytrop

! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
            laysolfr = min(lay+1,laytrop)
         o2cont = 4.35e-4_rb*colo2(lay)/(350.0_rb*2.0_rb)
         speccomb = colh2o(lay) + o2adj*strrat*colo2(lay)
         specparm = colh2o(lay)/speccomb 
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
!         odadj = specparm + o2adj * (1._rb - specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(22) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(22) + js
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl
         do ig = 1, ng22
            taug_rrtmg(lay,ngs23+ig) = speccomb * &
                (fac000 * absa(ind0,ig) + &
                 fac100 * absa(ind0+1,ig) + &
                 fac010 * absa(ind0+9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1,ig) + &
                 fac101 * absa(ind1+1,ig) + &
                 fac011 * absa(ind1+9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
                 colh2o(lay) * &
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                  (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) &
                 + o2cont
!            ssa(lay,ngs23+ig) = tauray/taug_rrtmg(lay,ngs23+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs23+ig) = sfluxref(ig,js) &
!                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taur_rrtmg(lay,ngs23+ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         o2cont = 4.35e-4_rb*colo2(lay)/(350.0_rb*2.0_rb)
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(22) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(22) + 1
!         tauray = colmol(lay) * rayl
         do ig = 1, ng22
            taug_rrtmg(lay,ngs23+ig) = colo2(lay) * o2adj * &
                (fac00(lay) * absb(ind0,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1,ig) + &
                 fac11(lay) * absb(ind1+1,ig)) + &
                 o2cont
!            ssa(lay,ngs23+ig) = tauray/taug_rrtmg(lay,ngs23+ig)
!            taur_rrtmg(lay,ngs23+ig) = tauray
         enddo
      enddo
      kout=kout + ng22
      end subroutine taumol22

!----------------------------------------------------------------------------
      subroutine taumol23(kout)
!----------------------------------------------------------------------------
!
!     band 23:  8050-12850 cm-1 (low - h2o; high - nothing)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng23, ngs24m
      use rrsw_kg23, only : absa, ka, forref, selfref, &
                            sfluxref, rayl

! ------- Declarations -------
      integer(kind=im), intent(inout):: kout

! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray, givfac

! Average Giver et al. correction factor for this band.
      givfac = 1.029_rb

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      layreffr = 6
      laysolfr = laytrop

! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
            laysolfr = min(lay+1,laytrop)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(23) + 1
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(23) + 1
         inds = indself(lay)
         indf = indfor(lay)
         do ig = 1, ng23
 !           tauray = colmol(lay) * rayl(ig)
            taug_rrtmg(lay,ngs24m+ig) = colh2o(lay) * &
                (givfac * (fac00(lay) * absa(ind0,ig) + &
                 fac10(lay) * absa(ind0+1,ig) + &
                 fac01(lay) * absa(ind1,ig) + &
                 fac11(lay) * absa(ind1+1,ig)) + &
                 selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + &
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) 
!            ssa(lay,ngs22+ig) = tauray/taug_rrtmg(lay,ngs22+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs24m+ig) = sfluxref(ig) 
!            taur_rrtmg(lay,ngs24m+ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         do ig = 1, ng23
!            taur_rrtmg(lay,ngs24m+ig) = colmol(lay) * rayl(ig)
!            ssa(lay,ngs22+ig) = 1.0_rb
            taug_rrtmg(lay, ngs24m+ig) = 0._rb
         enddo
      enddo
      kout= kout + ng23
      end subroutine taumol23

!----------------------------------------------------------------------------
      subroutine taumol24(taux,sfluxx)
!----------------------------------------------------------------------------
!
!     band 24:  12850-16000 cm-1 (low - h2o, o2; high - o2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng24
      use rrsw_kg24, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, abso3a, abso3b, rayla, raylb

! ------- Declarations -------
      real(kind=rb), intent(out) ::  taux(nlayers, ng24), sfluxx(ng24)
!taurx(nlayers, ng24), 
! Local
      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, layreffr
      real (kind=rb)  ::  fac000, fac001, fac010, fac011, fac100, fac101
      real (kind=rb)  ::  fac110, fac111, fs, speccomb, specmult, specparm
      real (kind=rb)  ::  tauray, strrat
! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.
      taux(:,:)=0.d0  
      strrat = 0.124692_rb
      layreffr = 1
      laysolfr = laytrop
! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
            laysolfr = min(lay+1, laytrop)
         speccomb = colh2o(lay) + strrat*colo2(lay)
         specparm = colh2o(lay)/speccomb
         if (specparm .ge. oneminus) specparm = oneminus
         specmult = 8._rb*(specparm)
         js = 1 + int(specmult)
         fs = mod(specmult, 1._rb )
         fac000 = (1._rb - fs) * fac00(lay)
         fac010 = (1._rb - fs) * fac10(lay)
         fac100 = fs * fac00(lay)
         fac110 = fs * fac10(lay)
         fac001 = (1._rb - fs) * fac01(lay)
         fac011 = (1._rb - fs) * fac11(lay)
         fac101 = fs * fac01(lay)
         fac111 = fs * fac11(lay)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(24) + js
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(24) + js
         inds = indself(lay)
         indf = indfor(lay)
         do ig = 1, ng24
            taux(lay,ig) = speccomb * &
                (fac000 * absa(ind0,ig) + &
                 fac100 * absa(ind0+1,ig) + &
                 fac010 * absa(ind0+9,ig) + &
                 fac110 * absa(ind0+10,ig) + &
                 fac001 * absa(ind1,ig) + &
                 fac101 * absa(ind1+1,ig) + &
                 fac011 * absa(ind1+9,ig) + &
                 fac111 * absa(ind1+10,ig)) + &
!                 colo3(lay) * abso3a(ig) + & 
                 colh2o(lay) * & 
                 (selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + & 
                 forfrac(lay) * &
                (forref(indf+1,ig) - forref(indf,ig))))
!            ssa(lay,ig) = tauray/taug_rrtmg(lay,ngs23+ig)
            if (lay .eq. laysolfr) sfluxx(ig) = sfluxref(ig,js) &
                 + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
!            taurx(lay, ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(24) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(24) + 1
         do ig = 1, ng24
!            tauray = colmol(lay) * raylb(ig)
            taux(lay,ig) = colo2(lay) * &
                (fac00(lay) * absb(ind0,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1,ig) + &
                 fac11(lay) * absb(ind1+1,ig)) !+ &
!                 colo3(lay) * abso3b(ig)
! 
!            ssa(lay,ngs23+ig) = tauray/taug_rrtmg(lay,ngs23+ig)
!            taurx(lay,ig) = tauray
         enddo
      enddo


      end subroutine taumol24

!----------------------------------------------------------------------------
      subroutine taumol25(taux, sfluxx)
!----------------------------------------------------------------------------
!
!     band 25:  16000-22650 cm-1 (low - h2o; high - nothing)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng25
      use rrsw_kg25, only : absa, ka, &
                            sfluxref, abso3a, abso3b, rayl

! ------- Declarations -------

      real(kind=rb), intent(out) ::  taux(nlayers,ng25), sfluxx(ng25)
!, taurx(nlayers, ng25),&  
! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray

! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  
      taux(:,:)=0.d0
      layreffr = 2
      laysolfr = laytrop      
! Lower atmosphere loop
      do lay = 1, laytrop
         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
              laysolfr = min(lay+1,laytrop)
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(25) + 1
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(25) + 1
         do ig = 1, ng25
!           tauray = colmol(lay) * rayl(ig)
            taux(lay,ig) = colh2o(lay) * &
                 (fac00(lay) * absa(ind0,ig) + &
                 fac10(lay) * absa(ind0+1,ig) + &
                 fac01(lay) * absa(ind1,ig) + &
                 fac11(lay) * absa(ind1+1,ig))!+&
!                 colo3(lay) * abso3a(ig) 
            if (lay .eq. laysolfr) sfluxx(ig) = sfluxref(ig) 
!            taurx(lay,ig) = tauray
         enddo
      enddo
! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         do ig = 1, ng25
!            tauray = colmol(lay) * rayl(ig)
             taux(lay,ig)= 1.d-30
!            taux(lay,ig) = colo3(lay) * abso3b(ig) 
!            ssa(lay,ngs24+ig) = tauray/taug_rrtmg(lay,ngs24+ig)
!            taurx(lay,ig) = tauray
         enddo
      enddo

!      do lay= nlayers,1,-1
!      write(110,'(10X, 2E10.3)') colo3(lay), colh2o(lay)
!      enddo

!      write(112,'(6F10.3)')(sfluxx(ig),ig=1, ng25)
!      do lay=nlayers,1,-1
!      write(112,'(6E10.3)')(taux(lay,ig),ig=1, ng25)
!      enddo

      end subroutine taumol25

!----------------------------------------------------------------------------
      subroutine taumol29(kout)
!----------------------------------------------------------------------------
!
!     band 29:  820-2600 cm-1 (low - h2o; high - co2)
!
!----------------------------------------------------------------------------

! ------- Modules -------

      use parrrsw, only : ng29, ngs16
      use rrsw_kg29, only : absa, ka, absb, kb, forref, selfref, &
                            sfluxref, absh2o, absco2, rayl
! ------- Declarations -------
     integer(kind=im),intent(inout):: kout
! Local

      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
                          layreffr
      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
                       fac110, fac111, fs, speccomb, specmult, specparm, &
                       tauray
! Compute the optical depth by interpolating in ln(pressure), 
! temperature, and appropriate species.  Below LAYTROP, the water
! vapor self-continuum is interpolated (in temperature) separately.  

      layreffr = 49

! Lower atmosphere loop
      do lay = 1, laytrop
         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(29) + 1
         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(29) + 1
         inds = indself(lay)
         indf = indfor(lay)
!         tauray = colmol(lay) * rayl

         do ig = 1, ng29
            taug_rrtmg(lay,ngs16+ig) = colh2o(lay) * &
               ((fac00(lay) * absa(ind0,ig) + &
                 fac10(lay) * absa(ind0+1,ig) + &
                 fac01(lay) * absa(ind1,ig) + &
                 fac11(lay) * absa(ind1+1,ig)) + &
                 selffac(lay) * (selfref(inds,ig) + &
                 selffrac(lay) * &
                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
                 forfac(lay) * (forref(indf,ig) + & 
                 forfrac(lay) * &
                 (forref(indf+1,ig) - forref(indf,ig)))) &
                 + colco2(lay) * absco2(ig) 
!            ssa(lay,ngs16+ig) = tauray/taug_rrtmg(lay,ngs16+ig)
!            taur_rrtmg(lay,ngs16+ig) = tauray
         enddo
      enddo

      laysolfr = nlayers

! Upper atmosphere loop
      do lay = laytrop+1, nlayers
         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
            laysolfr = lay
         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(29) + 1
         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(29) + 1
!         tauray = colmol(lay) * rayl

         do ig = 1, ng29
            taug_rrtmg(lay,ngs16+ig) = colco2(lay) * &
                (fac00(lay) * absb(ind0,ig) + &
                 fac10(lay) * absb(ind0+1,ig) + &
                 fac01(lay) * absb(ind1,ig) + &
                 fac11(lay) * absb(ind1+1,ig)) &  
                 + colh2o(lay) * absh2o(ig) 
!            ssa(lay,ngs16+ig) = tauray/taug_rrtmg(lay,ngs16+ig)
!            if (lay .eq. laysolfr) sfluxzen(ngs16+ig) = sfluxref(ig) 
!            taur_rrtmg(lay,ngs16+ig) = tauray
         enddo
      enddo
      kout= kout + ng29  
      end subroutine taumol29

      end subroutine taumol_sw

      end module rrtmg_sw_taumol

