source: palm/trunk/LIB/rrtmg/rrtmg_lw_taumol.f90 @ 3812

Last change on this file since 3812 was 1585, checked in by maronga, 10 years ago

Added support for RRTMG radiation code

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