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