source: palm/trunk/LIB/rrtmg/rrtmg_sw_taumol.f90 @ 4210

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

Added support for RRTMG radiation code

File size: 61.9 KB
Line 
1!     path:      $Source$
2!     author:    $Author: mike $
3!     revision:  $Revision: 11661 $
4!     created:   $Date: 2009-05-22 18:22:22 -0400 (Fri, 22 May 2009) $
5
6      module rrtmg_sw_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 parrrsw, only : mg, jpband, nbndsw, ngptsw
22      use rrsw_con, only: oneminus
23      use rrsw_wvn, only: nspa, nspb
24      use rrsw_vsn, only: hvrtau, hnamtau
25
26      implicit none
27
28      contains
29
30!----------------------------------------------------------------------------
31      subroutine taumol_sw(nlayers, &
32                           colh2o, colco2, colch4, colo2, colo3, colmol, &
33                           laytrop, jp, jt, jt1, &
34                           fac00, fac01, fac10, fac11, &
35                           selffac, selffrac, indself, forfac, forfrac, indfor, &
36                           sfluxzen, taug, taur)
37!----------------------------------------------------------------------------
38
39! ******************************************************************************
40! *                                                                            *
41! *                 Optical depths developed for the                           *
42! *                                                                            *
43! *               RAPID RADIATIVE TRANSFER MODEL (RRTM)                        *
44! *                                                                            *
45! *                                                                            *
46! *           ATMOSPHERIC AND ENVIRONMENTAL RESEARCH, INC.                     *
47! *                       131 HARTWELL AVENUE                                  *
48! *                       LEXINGTON, MA 02421                                  *
49! *                                                                            *
50! *                                                                            *
51! *                          ELI J. MLAWER                                     *
52! *                        JENNIFER DELAMERE                                   *
53! *                        STEVEN J. TAUBMAN                                   *
54! *                        SHEPARD A. CLOUGH                                   *
55! *                                                                            *
56! *                                                                            *
57! *                                                                            *
58! *                                                                            *
59! *                      email:  mlawer@aer.com                                *
60! *                      email:  jdelamer@aer.com                              *
61! *                                                                            *
62! *       The authors wish to acknowledge the contributions of the             *
63! *       following people:  Patrick D. Brown, Michael J. Iacono,              *
64! *       Ronald E. Farren, Luke Chen, Robert Bergstrom.                       *
65! *                                                                            *
66! ******************************************************************************
67! *    TAUMOL                                                                  *
68! *                                                                            *
69! *    This file contains the subroutines TAUGBn (where n goes from            *
70! *    1 to 28).  TAUGBn calculates the optical depths and Planck fractions    *
71! *    per g-value and layer for band n.                                       *
72! *                                                                            *
73! * Output:  optical depths (unitless)                                         *
74! *          fractions needed to compute Planck functions at every layer       *
75! *              and g-value                                                   *
76! *                                                                            *
77! *    COMMON /TAUGCOM/  TAUG(MXLAY,MG)                                        *
78! *    COMMON /PLANKG/   FRACS(MXLAY,MG)                                       *
79! *                                                                            *
80! * Input                                                                      *
81! *                                                                            *
82! *    PARAMETER (MG=16, MXLAY=203, NBANDS=14)                                 *
83! *                                                                            *
84! *    COMMON /FEATURES/ NG(NBANDS),NSPA(NBANDS),NSPB(NBANDS)                  *
85! *    COMMON /PRECISE/  ONEMINUS                                              *
86! *    COMMON /PROFILE/  NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),                    *
87! *   &                  PZ(0:MXLAY),TZ(0:MXLAY),TBOUND                        *
88! *    COMMON /PROFDATA/ LAYTROP,LAYSWTCH,LAYLOW,                              *
89! *   &                  COLH2O(MXLAY),COLCO2(MXLAY),                          *
90! *   &                  COLO3(MXLAY),COLN2O(MXLAY),COLCH4(MXLAY),             *
91! *   &                  COLO2(MXLAY),CO2MULT(MXLAY)                           *
92! *    COMMON /INTFAC/   FAC00(MXLAY),FAC01(MXLAY),                            *
93! *   &                  FAC10(MXLAY),FAC11(MXLAY)                             *
94! *    COMMON /INTIND/   JP(MXLAY),JT(MXLAY),JT1(MXLAY)                        *
95! *    COMMON /SELF/     SELFFAC(MXLAY), SELFFRAC(MXLAY), INDSELF(MXLAY)       *
96! *                                                                            *
97! *    Description:                                                            *
98! *    NG(IBAND) - number of g-values in band IBAND                            *
99! *    NSPA(IBAND) - for the lower atmosphere, the number of reference         *
100! *                  atmospheres that are stored for band IBAND per            *
101! *                  pressure level and temperature.  Each of these            *
102! *                  atmospheres has different relative amounts of the         *
103! *                  key species for the band (i.e. different binary           *
104! *                  species parameters).                                      *
105! *    NSPB(IBAND) - same for upper atmosphere                                 *
106! *    ONEMINUS - since problems are caused in some cases by interpolation     *
107! *               parameters equal to or greater than 1, for these cases       *
108! *               these parameters are set to this value, slightly < 1.        *
109! *    PAVEL - layer pressures (mb)                                            *
110! *    TAVEL - layer temperatures (degrees K)                                  *
111! *    PZ - level pressures (mb)                                               *
112! *    TZ - level temperatures (degrees K)                                     *
113! *    LAYTROP - layer at which switch is made from one combination of         *
114! *              key species to another                                        *
115! *    COLH2O, COLCO2, COLO3, COLN2O, COLCH4 - column amounts of water         *
116! *              vapor,carbon dioxide, ozone, nitrous ozide, methane,          *
117! *              respectively (molecules/cm**2)                                *
118! *    CO2MULT - for bands in which carbon dioxide is implemented as a         *
119! *              trace species, this is the factor used to multiply the        *
120! *              band's average CO2 absorption coefficient to get the added    *
121! *              contribution to the optical depth relative to 355 ppm.        *
122! *    FACij(LAY) - for layer LAY, these are factors that are needed to        *
123! *                 compute the interpolation factors that multiply the        *
124! *                 appropriate reference k-values.  A value of 0 (1) for      *
125! *                 i,j indicates that the corresponding factor multiplies     *
126! *                 reference k-value for the lower (higher) of the two        *
127! *                 appropriate temperatures, and altitudes, respectively.     *
128! *    JP - the index of the lower (in altitude) of the two appropriate        *
129! *         reference pressure levels needed for interpolation                 *
130! *    JT, JT1 - the indices of the lower of the two appropriate reference     *
131! *              temperatures needed for interpolation (for pressure           *
132! *              levels JP and JP+1, respectively)                             *
133! *    SELFFAC - scale factor needed to water vapor self-continuum, equals     *
134! *              (water vapor density)/(atmospheric density at 296K and        *
135! *              1013 mb)                                                      *
136! *    SELFFRAC - factor needed for temperature interpolation of reference     *
137! *               water vapor self-continuum data                              *
138! *    INDSELF - index of the lower of the two appropriate reference           *
139! *              temperatures needed for the self-continuum interpolation      *
140! *                                                                            *
141! * Data input                                                                 *
142! *    COMMON /Kn/ KA(NSPA(n),5,13,MG), KB(NSPB(n),5,13:59,MG), SELFREF(10,MG) *
143! *       (note:  n is the band number)                                        *
144! *                                                                            *
145! *    Description:                                                            *
146! *    KA - k-values for low reference atmospheres (no water vapor             *
147! *         self-continuum) (units: cm**2/molecule)                            *
148! *    KB - k-values for high reference atmospheres (all sources)              *
149! *         (units: cm**2/molecule)                                            *
150! *    SELFREF - k-values for water vapor self-continuum for reference         *
151! *              atmospheres (used below LAYTROP)                              *
152! *              (units: cm**2/molecule)                                       *
153! *                                                                            *
154! *    DIMENSION ABSA(65*NSPA(n),MG), ABSB(235*NSPB(n),MG)                     *
155! *    EQUIVALENCE (KA,ABSA),(KB,ABSB)                                         *
156! *                                                                            *
157! *****************************************************************************
158!
159! Modifications
160!
161! Revised: Adapted to F90 coding, J.-J.Morcrette, ECMWF, Feb 2003
162! Revised: Modified for g-point reduction, MJIacono, AER, Dec 2003
163! Revised: Reformatted for consistency with rrtmg_lw, MJIacono, AER, Jul 2006
164!
165! ------- Declarations -------
166
167! ----- Input -----
168      integer(kind=im), intent(in) :: nlayers            ! total number of layers
169
170      integer(kind=im), intent(in) :: laytrop            ! tropopause layer index
171      integer(kind=im), intent(in) :: jp(:)              !
172                                                         !   Dimensions: (nlayers)
173      integer(kind=im), intent(in) :: jt(:)              !
174                                                         !   Dimensions: (nlayers)
175      integer(kind=im), intent(in) :: jt1(:)             !
176                                                         !   Dimensions: (nlayers)
177
178      real(kind=rb), intent(in) :: colh2o(:)             ! column amount (h2o)
179                                                         !   Dimensions: (nlayers)
180      real(kind=rb), intent(in) :: colco2(:)             ! column amount (co2)
181                                                         !   Dimensions: (nlayers)
182      real(kind=rb), intent(in) :: colo3(:)              ! column amount (o3)
183                                                         !   Dimensions: (nlayers)
184      real(kind=rb), intent(in) :: colch4(:)             ! column amount (ch4)
185                                                         !   Dimensions: (nlayers)
186                                                         !   Dimensions: (nlayers)
187      real(kind=rb), intent(in) :: colo2(:)              ! column amount (o2)
188                                                         !   Dimensions: (nlayers)
189      real(kind=rb), intent(in) :: colmol(:)             !
190                                                         !   Dimensions: (nlayers)
191
192      integer(kind=im), intent(in) :: indself(:)   
193                                                         !   Dimensions: (nlayers)
194      integer(kind=im), intent(in) :: indfor(:)
195                                                         !   Dimensions: (nlayers)
196      real(kind=rb), intent(in) :: selffac(:)
197                                                         !   Dimensions: (nlayers)
198      real(kind=rb), intent(in) :: selffrac(:)
199                                                         !   Dimensions: (nlayers)
200      real(kind=rb), intent(in) :: forfac(:)
201                                                         !   Dimensions: (nlayers)
202      real(kind=rb), intent(in) :: forfrac(:)
203                                                         !   Dimensions: (nlayers)
204
205      real(kind=rb), intent(in) :: &                     !
206                       fac00(:), fac01(:), &             !   Dimensions: (nlayers)
207                       fac10(:), fac11(:) 
208
209! ----- Output -----
210      real(kind=rb), intent(out) :: sfluxzen(:)          ! solar source function
211                                                         !   Dimensions: (ngptsw)
212      real(kind=rb), intent(out) :: taug(:,:)            ! gaseous optical depth
213                                                         !   Dimensions: (nlayers,ngptsw)
214      real(kind=rb), intent(out) :: taur(:,:)            ! Rayleigh
215                                                         !   Dimensions: (nlayers,ngptsw)
216!      real(kind=rb), intent(out) :: ssa(:,:)            ! single scattering albedo (inactive)
217                                                         !   Dimensions: (nlayers,ngptsw)
218
219      hvrtau = '$Revision: 11661 $'
220
221! Calculate gaseous optical depth and planck fractions for each spectral band.
222
223      call taumol16
224      call taumol17
225      call taumol18
226      call taumol19
227      call taumol20
228      call taumol21
229      call taumol22
230      call taumol23
231      call taumol24
232      call taumol25
233      call taumol26
234      call taumol27
235      call taumol28
236      call taumol29
237
238!-------------
239      contains
240!-------------
241
242!----------------------------------------------------------------------------
243      subroutine taumol16
244!----------------------------------------------------------------------------
245!
246!     band 16:  2600-3250 cm-1 (low - h2o,ch4; high - ch4)
247!
248!----------------------------------------------------------------------------
249
250! ------- Modules -------
251
252      use parrrsw, only : ng16
253      use rrsw_kg16, only : absa, ka, absb, kb, forref, selfref, &
254                            sfluxref, rayl
255
256! ------- Declarations -------
257
258! Local
259
260      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
261                          layreffr
262      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
263                       fac110, fac111, fs, speccomb, specmult, specparm, &
264                       tauray, strrat1
265
266! Compute the optical depth by interpolating in ln(pressure),
267! temperature, and appropriate species.  Below LAYTROP, the water
268! vapor self-continuum is interpolated (in temperature) separately. 
269
270      strrat1 = 252.131_rb
271      layreffr = 18
272
273! Lower atmosphere loop
274      do lay = 1, laytrop
275         speccomb = colh2o(lay) + strrat1*colch4(lay)
276         specparm = colh2o(lay)/speccomb 
277         if (specparm .ge. oneminus) specparm = oneminus
278         specmult = 8._rb*(specparm)
279         js = 1 + int(specmult)
280         fs = mod(specmult, 1._rb )
281         fac000 = (1._rb - fs) * fac00(lay)
282         fac010 = (1._rb - fs) * fac10(lay)
283         fac100 = fs * fac00(lay)
284         fac110 = fs * fac10(lay)
285         fac001 = (1._rb - fs) * fac01(lay)
286         fac011 = (1._rb - fs) * fac11(lay)
287         fac101 = fs * fac01(lay)
288         fac111 = fs * fac11(lay)
289         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(16) + js
290         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(16) + js
291         inds = indself(lay)
292         indf = indfor(lay)
293         tauray = colmol(lay) * rayl
294
295         do ig = 1, ng16
296            taug(lay,ig) = speccomb * &
297                (fac000 * absa(ind0   ,ig) + &
298                 fac100 * absa(ind0 +1,ig) + &
299                 fac010 * absa(ind0 +9,ig) + &
300                 fac110 * absa(ind0+10,ig) + &
301                 fac001 * absa(ind1   ,ig) + &
302                 fac101 * absa(ind1 +1,ig) + &
303                 fac011 * absa(ind1 +9,ig) + &
304                 fac111 * absa(ind1+10,ig)) + &
305                 colh2o(lay) * &
306                 (selffac(lay) * (selfref(inds,ig) + &
307                 selffrac(lay) * &
308                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
309                 forfac(lay) * (forref(indf,ig) + &
310                 forfrac(lay) * &
311                 (forref(indf+1,ig) - forref(indf,ig)))) 
312!            ssa(lay,ig) = tauray/taug(lay,ig)
313            taur(lay,ig) = tauray
314         enddo
315      enddo
316
317      laysolfr = nlayers
318
319! Upper atmosphere loop
320      do lay = laytrop+1, nlayers
321         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
322            laysolfr = lay
323         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(16) + 1
324         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(16) + 1
325         tauray = colmol(lay) * rayl
326
327         do ig = 1, ng16
328            taug(lay,ig) = colch4(lay) * &
329                (fac00(lay) * absb(ind0  ,ig) + &
330                 fac10(lay) * absb(ind0+1,ig) + &
331                 fac01(lay) * absb(ind1  ,ig) + &
332                 fac11(lay) * absb(ind1+1,ig)) 
333!            ssa(lay,ig) = tauray/taug(lay,ig)
334            if (lay .eq. laysolfr) sfluxzen(ig) = sfluxref(ig) 
335            taur(lay,ig) = tauray 
336         enddo
337      enddo
338
339      end subroutine taumol16
340
341!----------------------------------------------------------------------------
342      subroutine taumol17
343!----------------------------------------------------------------------------
344!
345!     band 17:  3250-4000 cm-1 (low - h2o,co2; high - h2o,co2)
346!
347!----------------------------------------------------------------------------
348
349! ------- Modules -------
350
351      use parrrsw, only : ng17, ngs16
352      use rrsw_kg17, only : absa, ka, absb, kb, forref, selfref, &
353                            sfluxref, rayl
354
355! ------- Declarations -------
356
357! Local
358
359      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
360                          layreffr
361      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
362                       fac110, fac111, fs, speccomb, specmult, specparm, &
363                       tauray, strrat
364
365! Compute the optical depth by interpolating in ln(pressure),
366! temperature, and appropriate species.  Below LAYTROP, the water
367! vapor self-continuum is interpolated (in temperature) separately. 
368
369      strrat = 0.364641_rb
370      layreffr = 30
371
372! Lower atmosphere loop
373      do lay = 1, laytrop
374         speccomb = colh2o(lay) + strrat*colco2(lay)
375         specparm = colh2o(lay)/speccomb 
376         if (specparm .ge. oneminus) specparm = oneminus
377         specmult = 8._rb*(specparm)
378         js = 1 + int(specmult)
379         fs = mod(specmult, 1._rb )
380         fac000 = (1._rb - fs) * fac00(lay)
381         fac010 = (1._rb - fs) * fac10(lay)
382         fac100 = fs * fac00(lay)
383         fac110 = fs * fac10(lay)
384         fac001 = (1._rb - fs) * fac01(lay)
385         fac011 = (1._rb - fs) * fac11(lay)
386         fac101 = fs * fac01(lay)
387         fac111 = fs * fac11(lay)
388         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(17) + js
389         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(17) + js
390         inds = indself(lay)
391         indf = indfor(lay)
392         tauray = colmol(lay) * rayl
393
394         do ig = 1, ng17
395            taug(lay,ngs16+ig) = speccomb * &
396                (fac000 * absa(ind0,ig) + &
397                 fac100 * absa(ind0+1,ig) + &
398                 fac010 * absa(ind0+9,ig) + &
399                 fac110 * absa(ind0+10,ig) + &
400                 fac001 * absa(ind1,ig) + &
401                 fac101 * absa(ind1+1,ig) + &
402                 fac011 * absa(ind1+9,ig) + &
403                 fac111 * absa(ind1+10,ig)) + &
404                 colh2o(lay) * &
405                 (selffac(lay) * (selfref(inds,ig) + &
406                 selffrac(lay) * &
407                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
408                 forfac(lay) * (forref(indf,ig) + &
409                 forfrac(lay) * &
410                 (forref(indf+1,ig) - forref(indf,ig)))) 
411!             ssa(lay,ngs16+ig) = tauray/taug(lay,ngs16+ig)
412            taur(lay,ngs16+ig) = tauray
413         enddo
414      enddo
415
416      laysolfr = nlayers
417
418! Upper atmosphere loop
419      do lay = laytrop+1, nlayers
420         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
421            laysolfr = lay
422         speccomb = colh2o(lay) + strrat*colco2(lay)
423         specparm = colh2o(lay)/speccomb 
424         if (specparm .ge. oneminus) specparm = oneminus
425         specmult = 4._rb*(specparm)
426         js = 1 + int(specmult)
427         fs = mod(specmult, 1._rb )
428         fac000 = (1._rb - fs) * fac00(lay)
429         fac010 = (1._rb - fs) * fac10(lay)
430         fac100 = fs * fac00(lay)
431         fac110 = fs * fac10(lay)
432         fac001 = (1._rb - fs) * fac01(lay)
433         fac011 = (1._rb - fs) * fac11(lay)
434         fac101 = fs * fac01(lay)
435         fac111 = fs * fac11(lay)
436         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(17) + js
437         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(17) + js
438         indf = indfor(lay)
439         tauray = colmol(lay) * rayl
440
441         do ig = 1, ng17
442            taug(lay,ngs16+ig) = speccomb * &
443                (fac000 * absb(ind0,ig) + &
444                 fac100 * absb(ind0+1,ig) + &
445                 fac010 * absb(ind0+5,ig) + &
446                 fac110 * absb(ind0+6,ig) + &
447                 fac001 * absb(ind1,ig) + &
448                 fac101 * absb(ind1+1,ig) + &
449                 fac011 * absb(ind1+5,ig) + &
450                 fac111 * absb(ind1+6,ig)) + &
451                 colh2o(lay) * &
452                 forfac(lay) * (forref(indf,ig) + &
453                 forfrac(lay) * &
454                 (forref(indf+1,ig) - forref(indf,ig))) 
455!            ssa(lay,ngs16+ig) = tauray/taug(lay,ngs16+ig)
456            if (lay .eq. laysolfr) sfluxzen(ngs16+ig) = sfluxref(ig,js) &
457               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
458            taur(lay,ngs16+ig) = tauray
459         enddo
460      enddo
461
462      end subroutine taumol17
463
464!----------------------------------------------------------------------------
465      subroutine taumol18
466!----------------------------------------------------------------------------
467!
468!     band 18:  4000-4650 cm-1 (low - h2o,ch4; high - ch4)
469!
470!----------------------------------------------------------------------------
471
472! ------- Modules -------
473
474      use parrrsw, only : ng18, ngs17
475      use rrsw_kg18, only : absa, ka, absb, kb, forref, selfref, &
476                            sfluxref, rayl
477
478! ------- Declarations -------
479
480! Local
481
482      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
483                          layreffr
484      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
485                       fac110, fac111, fs, speccomb, specmult, specparm, &
486                       tauray, strrat
487
488! Compute the optical depth by interpolating in ln(pressure),
489! temperature, and appropriate species.  Below LAYTROP, the water
490! vapor self-continuum is interpolated (in temperature) separately. 
491
492      strrat = 38.9589_rb
493      layreffr = 6
494      laysolfr = laytrop
495     
496! Lower atmosphere loop
497      do lay = 1, laytrop
498         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
499            laysolfr = min(lay+1,laytrop)
500         speccomb = colh2o(lay) + strrat*colch4(lay)
501         specparm = colh2o(lay)/speccomb 
502         if (specparm .ge. oneminus) specparm = oneminus
503         specmult = 8._rb*(specparm)
504         js = 1 + int(specmult)
505         fs = mod(specmult, 1._rb )
506         fac000 = (1._rb - fs) * fac00(lay)
507         fac010 = (1._rb - fs) * fac10(lay)
508         fac100 = fs * fac00(lay)
509         fac110 = fs * fac10(lay)
510         fac001 = (1._rb - fs) * fac01(lay)
511         fac011 = (1._rb - fs) * fac11(lay)
512         fac101 = fs * fac01(lay)
513         fac111 = fs * fac11(lay)
514         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(18) + js
515         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(18) + js
516         inds = indself(lay)
517         indf = indfor(lay)
518         tauray = colmol(lay) * rayl
519
520         do ig = 1, ng18
521            taug(lay,ngs17+ig) = speccomb * &
522                (fac000 * absa(ind0,ig) + &
523                 fac100 * absa(ind0+1,ig) + &
524                 fac010 * absa(ind0+9,ig) + &
525                 fac110 * absa(ind0+10,ig) + &
526                 fac001 * absa(ind1,ig) + &
527                 fac101 * absa(ind1+1,ig) + &
528                 fac011 * absa(ind1+9,ig) + &
529                 fac111 * absa(ind1+10,ig)) + &
530                 colh2o(lay) * &
531                 (selffac(lay) * (selfref(inds,ig) + &
532                 selffrac(lay) * &
533                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
534                 forfac(lay) * (forref(indf,ig) + &
535                 forfrac(lay) * &
536                 (forref(indf+1,ig) - forref(indf,ig)))) 
537!            ssa(lay,ngs17+ig) = tauray/taug(lay,ngs17+ig)
538            if (lay .eq. laysolfr) sfluxzen(ngs17+ig) = sfluxref(ig,js) &
539               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
540            taur(lay,ngs17+ig) = tauray
541         enddo
542      enddo
543
544! Upper atmosphere loop
545      do lay = laytrop+1, nlayers
546         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(18) + 1
547         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(18) + 1
548         tauray = colmol(lay) * rayl
549
550         do ig = 1, ng18
551            taug(lay,ngs17+ig) = colch4(lay) * &
552                (fac00(lay) * absb(ind0,ig) + &
553                 fac10(lay) * absb(ind0+1,ig) + &
554                 fac01(lay) * absb(ind1,ig) + &   
555                 fac11(lay) * absb(ind1+1,ig)) 
556!           ssa(lay,ngs17+ig) = tauray/taug(lay,ngs17+ig)
557           taur(lay,ngs17+ig) = tauray
558         enddo
559       enddo
560
561       end subroutine taumol18
562
563!----------------------------------------------------------------------------
564      subroutine taumol19
565!----------------------------------------------------------------------------
566!
567!     band 19:  4650-5150 cm-1 (low - h2o,co2; high - co2)
568!
569!----------------------------------------------------------------------------
570
571! ------- Modules -------
572
573      use parrrsw, only : ng19, ngs18
574      use rrsw_kg19, only : absa, ka, absb, kb, forref, selfref, &
575                            sfluxref, rayl
576
577! ------- Declarations -------
578
579! Local
580
581      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
582                          layreffr
583      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
584                       fac110, fac111, fs, speccomb, specmult, specparm, &
585                       tauray, strrat
586
587! Compute the optical depth by interpolating in ln(pressure),
588! temperature, and appropriate species.  Below LAYTROP, the water
589! vapor self-continuum is interpolated (in temperature) separately. 
590
591      strrat = 5.49281_rb
592      layreffr = 3
593      laysolfr = laytrop
594
595! Lower atmosphere loop     
596      do lay = 1, laytrop
597         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
598            laysolfr = min(lay+1,laytrop)
599         speccomb = colh2o(lay) + strrat*colco2(lay)
600         specparm = colh2o(lay)/speccomb 
601         if (specparm .ge. oneminus) specparm = oneminus
602         specmult = 8._rb*(specparm)
603         js = 1 + int(specmult)
604         fs = mod(specmult, 1._rb )
605         fac000 = (1._rb - fs) * fac00(lay)
606         fac010 = (1._rb - fs) * fac10(lay)
607         fac100 = fs * fac00(lay)
608         fac110 = fs * fac10(lay)
609         fac001 = (1._rb - fs) * fac01(lay)
610         fac011 = (1._rb - fs) * fac11(lay)
611         fac101 = fs * fac01(lay)
612         fac111 = fs * fac11(lay)
613         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(19) + js
614         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(19) + js
615         inds = indself(lay)
616         indf = indfor(lay)
617         tauray = colmol(lay) * rayl
618
619         do ig = 1 , ng19
620            taug(lay,ngs18+ig) = speccomb * &
621                (fac000 * absa(ind0,ig) + &
622                 fac100 * absa(ind0+1,ig) + &
623                 fac010 * absa(ind0+9,ig) + &
624                 fac110 * absa(ind0+10,ig) + &
625                 fac001 * absa(ind1,ig) + &
626                 fac101 * absa(ind1+1,ig) + &
627                 fac011 * absa(ind1+9,ig) + &
628                 fac111 * absa(ind1+10,ig)) + &
629                 colh2o(lay) * &
630                 (selffac(lay) * (selfref(inds,ig) + &
631                 selffrac(lay) * &
632                 (selfref(inds+1,ig) - selfref(inds,ig))) + & 
633                 forfac(lay) * (forref(indf,ig) + &
634                 forfrac(lay) * &
635                 (forref(indf+1,ig) - forref(indf,ig)))) 
636!            ssa(lay,ngs18+ig) = tauray/taug(lay,ngs18+ig)
637            if (lay .eq. laysolfr) sfluxzen(ngs18+ig) = sfluxref(ig,js) &
638               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
639            taur(lay,ngs18+ig) = tauray   
640         enddo
641      enddo
642
643! Upper atmosphere loop
644      do lay = laytrop+1, nlayers
645         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(19) + 1
646         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(19) + 1
647         tauray = colmol(lay) * rayl
648
649         do ig = 1 , ng19
650            taug(lay,ngs18+ig) = colco2(lay) * &
651                (fac00(lay) * absb(ind0,ig) + &
652                 fac10(lay) * absb(ind0+1,ig) + &
653                 fac01(lay) * absb(ind1,ig) + &
654                 fac11(lay) * absb(ind1+1,ig)) 
655!            ssa(lay,ngs18+ig) = tauray/taug(lay,ngs18+ig)
656            taur(lay,ngs18+ig) = tauray   
657         enddo
658      enddo
659
660      end subroutine taumol19
661
662!----------------------------------------------------------------------------
663      subroutine taumol20
664!----------------------------------------------------------------------------
665!
666!     band 20:  5150-6150 cm-1 (low - h2o; high - h2o)
667!
668!----------------------------------------------------------------------------
669
670! ------- Modules -------
671
672      use parrrsw, only : ng20, ngs19
673      use rrsw_kg20, only : absa, ka, absb, kb, forref, selfref, &
674                            sfluxref, absch4, rayl
675
676      implicit none
677
678! ------- Declarations -------
679
680! Local
681
682      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
683                          layreffr
684      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
685                       fac110, fac111, fs, speccomb, specmult, specparm, &
686                       tauray
687
688! Compute the optical depth by interpolating in ln(pressure),
689! temperature, and appropriate species.  Below LAYTROP, the water
690! vapor self-continuum is interpolated (in temperature) separately. 
691
692      layreffr = 3
693      laysolfr = laytrop
694
695! Lower atmosphere loop
696      do lay = 1, laytrop
697         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
698            laysolfr = min(lay+1,laytrop)
699         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(20) + 1
700         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(20) + 1
701         inds = indself(lay)
702         indf = indfor(lay)
703         tauray = colmol(lay) * rayl
704
705         do ig = 1, ng20
706            taug(lay,ngs19+ig) = colh2o(lay) * &
707               ((fac00(lay) * absa(ind0,ig) + &
708                 fac10(lay) * absa(ind0+1,ig) + &
709                 fac01(lay) * absa(ind1,ig) + &
710                 fac11(lay) * absa(ind1+1,ig)) + &
711                 selffac(lay) * (selfref(inds,ig) + & 
712                 selffrac(lay) * &
713                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
714                 forfac(lay) * (forref(indf,ig) + &
715                 forfrac(lay) * &
716                 (forref(indf+1,ig) - forref(indf,ig)))) &
717                 + colch4(lay) * absch4(ig)
718!            ssa(lay,ngs19+ig) = tauray/taug(lay,ngs19+ig)
719            taur(lay,ngs19+ig) = tauray 
720            if (lay .eq. laysolfr) sfluxzen(ngs19+ig) = sfluxref(ig) 
721         enddo
722      enddo
723
724! Upper atmosphere loop
725      do lay = laytrop+1, nlayers
726         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(20) + 1
727         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(20) + 1
728         indf = indfor(lay)
729         tauray = colmol(lay) * rayl
730
731         do ig = 1, ng20
732            taug(lay,ngs19+ig) = colh2o(lay) * &
733                (fac00(lay) * absb(ind0,ig) + &
734                 fac10(lay) * absb(ind0+1,ig) + &
735                 fac01(lay) * absb(ind1,ig) + &
736                 fac11(lay) * absb(ind1+1,ig) + &
737                 forfac(lay) * (forref(indf,ig) + &
738                 forfrac(lay) * &
739                 (forref(indf+1,ig) - forref(indf,ig)))) + &
740                 colch4(lay) * absch4(ig)
741!            ssa(lay,ngs19+ig) = tauray/taug(lay,ngs19+ig)
742            taur(lay,ngs19+ig) = tauray 
743         enddo
744      enddo
745
746      end subroutine taumol20
747
748!----------------------------------------------------------------------------
749      subroutine taumol21
750!----------------------------------------------------------------------------
751!
752!     band 21:  6150-7700 cm-1 (low - h2o,co2; high - h2o,co2)
753!
754!----------------------------------------------------------------------------
755
756! ------- Modules -------
757
758      use parrrsw, only : ng21, ngs20
759      use rrsw_kg21, only : absa, ka, absb, kb, forref, selfref, &
760                            sfluxref, rayl
761
762! ------- Declarations -------
763
764! Local
765
766      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
767                          layreffr
768      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
769                       fac110, fac111, fs, speccomb, specmult, specparm, &
770                       tauray, strrat
771
772! Compute the optical depth by interpolating in ln(pressure),
773! temperature, and appropriate species.  Below LAYTROP, the water
774! vapor self-continuum is interpolated (in temperature) separately. 
775
776      strrat = 0.0045321_rb
777      layreffr = 8
778      laysolfr = laytrop
779     
780! Lower atmosphere loop
781      do lay = 1, laytrop
782         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
783            laysolfr = min(lay+1,laytrop)
784         speccomb = colh2o(lay) + strrat*colco2(lay)
785         specparm = colh2o(lay)/speccomb 
786         if (specparm .ge. oneminus) specparm = oneminus
787         specmult = 8._rb*(specparm)
788         js = 1 + int(specmult)
789         fs = mod(specmult, 1._rb )
790         fac000 = (1._rb - fs) * fac00(lay)
791         fac010 = (1._rb - fs) * fac10(lay)
792         fac100 = fs * fac00(lay)
793         fac110 = fs * fac10(lay)
794         fac001 = (1._rb - fs) * fac01(lay)
795         fac011 = (1._rb - fs) * fac11(lay)
796         fac101 = fs * fac01(lay)
797         fac111 = fs * fac11(lay)
798         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(21) + js
799         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(21) + js
800         inds = indself(lay)
801         indf = indfor(lay)
802         tauray = colmol(lay) * rayl
803
804         do ig = 1, ng21
805            taug(lay,ngs20+ig) = speccomb * &
806                (fac000 * absa(ind0,ig) + &
807                 fac100 * absa(ind0+1,ig) + &
808                 fac010 * absa(ind0+9,ig) + &
809                 fac110 * absa(ind0+10,ig) + &
810                 fac001 * absa(ind1,ig) + &
811                 fac101 * absa(ind1+1,ig) + &
812                 fac011 * absa(ind1+9,ig) + &
813                 fac111 * absa(ind1+10,ig)) + &
814                 colh2o(lay) * &
815                 (selffac(lay) * (selfref(inds,ig) + &
816                 selffrac(lay) * &
817                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
818                 forfac(lay) * (forref(indf,ig) + &
819                 forfrac(lay) * &
820                 (forref(indf+1,ig) - forref(indf,ig))))
821!            ssa(lay,ngs20+ig) = tauray/taug(lay,ngs20+ig)
822            if (lay .eq. laysolfr) sfluxzen(ngs20+ig) = sfluxref(ig,js) &
823               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
824            taur(lay,ngs20+ig) = tauray
825         enddo
826      enddo
827
828! Upper atmosphere loop
829      do lay = laytrop+1, nlayers
830         speccomb = colh2o(lay) + strrat*colco2(lay)
831         specparm = colh2o(lay)/speccomb 
832         if (specparm .ge. oneminus) specparm = oneminus
833         specmult = 4._rb*(specparm)
834         js = 1 + int(specmult)
835         fs = mod(specmult, 1._rb )
836         fac000 = (1._rb - fs) * fac00(lay)
837         fac010 = (1._rb - fs) * fac10(lay)
838         fac100 = fs * fac00(lay)
839         fac110 = fs * fac10(lay)
840         fac001 = (1._rb - fs) * fac01(lay)
841         fac011 = (1._rb - fs) * fac11(lay)
842         fac101 = fs * fac01(lay)
843         fac111 = fs * fac11(lay)
844         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(21) + js
845         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(21) + js
846         indf = indfor(lay)
847         tauray = colmol(lay) * rayl
848
849         do ig = 1, ng21
850            taug(lay,ngs20+ig) = speccomb * &
851                (fac000 * absb(ind0,ig) + &
852                 fac100 * absb(ind0+1,ig) + &
853                 fac010 * absb(ind0+5,ig) + &
854                 fac110 * absb(ind0+6,ig) + &
855                 fac001 * absb(ind1,ig) + &
856                 fac101 * absb(ind1+1,ig) + &
857                 fac011 * absb(ind1+5,ig) + &
858                 fac111 * absb(ind1+6,ig)) + &
859                 colh2o(lay) * &
860                 forfac(lay) * (forref(indf,ig) + &
861                 forfrac(lay) * &
862                 (forref(indf+1,ig) - forref(indf,ig)))
863!            ssa(lay,ngs20+ig) = tauray/taug(lay,ngs20+ig)
864            taur(lay,ngs20+ig) = tauray
865         enddo
866      enddo
867
868      end subroutine taumol21
869
870!----------------------------------------------------------------------------
871      subroutine taumol22
872!----------------------------------------------------------------------------
873!
874!     band 22:  7700-8050 cm-1 (low - h2o,o2; high - o2)
875!
876!----------------------------------------------------------------------------
877
878! ------- Modules -------
879
880      use parrrsw, only : ng22, ngs21
881      use rrsw_kg22, only : absa, ka, absb, kb, forref, selfref, &
882                            sfluxref, rayl
883
884! ------- Declarations -------
885
886! Local
887
888      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
889                          layreffr
890      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
891                       fac110, fac111, fs, speccomb, specmult, specparm, &
892                       tauray, o2adj, o2cont, strrat
893
894! The following factor is the ratio of total O2 band intensity (lines
895! and Mate continuum) to O2 band intensity (line only).  It is needed
896! to adjust the optical depths since the k's include only lines.
897      o2adj = 1.6_rb
898     
899! Compute the optical depth by interpolating in ln(pressure),
900! temperature, and appropriate species.  Below LAYTROP, the water
901! vapor self-continuum is interpolated (in temperature) separately. 
902
903      strrat = 0.022708_rb
904      layreffr = 2
905      laysolfr = laytrop
906
907! Lower atmosphere loop
908      do lay = 1, laytrop
909         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
910            laysolfr = min(lay+1,laytrop)
911         o2cont = 4.35e-4_rb*colo2(lay)/(350.0_rb*2.0_rb)
912         speccomb = colh2o(lay) + o2adj*strrat*colo2(lay)
913         specparm = colh2o(lay)/speccomb 
914         if (specparm .ge. oneminus) specparm = oneminus
915         specmult = 8._rb*(specparm)
916!         odadj = specparm + o2adj * (1._rb - specparm)
917         js = 1 + int(specmult)
918         fs = mod(specmult, 1._rb )
919         fac000 = (1._rb - fs) * fac00(lay)
920         fac010 = (1._rb - fs) * fac10(lay)
921         fac100 = fs * fac00(lay)
922         fac110 = fs * fac10(lay)
923         fac001 = (1._rb - fs) * fac01(lay)
924         fac011 = (1._rb - fs) * fac11(lay)
925         fac101 = fs * fac01(lay)
926         fac111 = fs * fac11(lay)
927         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(22) + js
928         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(22) + js
929         inds = indself(lay)
930         indf = indfor(lay)
931         tauray = colmol(lay) * rayl
932
933         do ig = 1, ng22
934            taug(lay,ngs21+ig) = speccomb * &
935                (fac000 * absa(ind0,ig) + &
936                 fac100 * absa(ind0+1,ig) + &
937                 fac010 * absa(ind0+9,ig) + &
938                 fac110 * absa(ind0+10,ig) + &
939                 fac001 * absa(ind1,ig) + &
940                 fac101 * absa(ind1+1,ig) + &
941                 fac011 * absa(ind1+9,ig) + &
942                 fac111 * absa(ind1+10,ig)) + &
943                 colh2o(lay) * &
944                 (selffac(lay) * (selfref(inds,ig) + &
945                 selffrac(lay) * &
946                  (selfref(inds+1,ig) - selfref(inds,ig))) + &
947                 forfac(lay) * (forref(indf,ig) + &
948                 forfrac(lay) * &
949                 (forref(indf+1,ig) - forref(indf,ig)))) &
950                 + o2cont
951!            ssa(lay,ngs21+ig) = tauray/taug(lay,ngs21+ig)
952            if (lay .eq. laysolfr) sfluxzen(ngs21+ig) = sfluxref(ig,js) &
953                + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
954            taur(lay,ngs21+ig) = tauray
955         enddo
956      enddo
957
958! Upper atmosphere loop
959      do lay = laytrop+1, nlayers
960         o2cont = 4.35e-4_rb*colo2(lay)/(350.0_rb*2.0_rb)
961         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(22) + 1
962         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(22) + 1
963         tauray = colmol(lay) * rayl
964
965         do ig = 1, ng22
966            taug(lay,ngs21+ig) = colo2(lay) * o2adj * &
967                (fac00(lay) * absb(ind0,ig) + &
968                 fac10(lay) * absb(ind0+1,ig) + &
969                 fac01(lay) * absb(ind1,ig) + &
970                 fac11(lay) * absb(ind1+1,ig)) + &
971                 o2cont
972!            ssa(lay,ngs21+ig) = tauray/taug(lay,ngs21+ig)
973            taur(lay,ngs21+ig) = tauray
974         enddo
975      enddo
976
977      end subroutine taumol22
978
979!----------------------------------------------------------------------------
980      subroutine taumol23
981!----------------------------------------------------------------------------
982!
983!     band 23:  8050-12850 cm-1 (low - h2o; high - nothing)
984!
985!----------------------------------------------------------------------------
986
987! ------- Modules -------
988
989      use parrrsw, only : ng23, ngs22
990      use rrsw_kg23, only : absa, ka, forref, selfref, &
991                            sfluxref, rayl
992
993! ------- Declarations -------
994
995! Local
996
997      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
998                          layreffr
999      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1000                       fac110, fac111, fs, speccomb, specmult, specparm, &
1001                       tauray, givfac
1002
1003! Average Giver et al. correction factor for this band.
1004      givfac = 1.029_rb
1005
1006! Compute the optical depth by interpolating in ln(pressure),
1007! temperature, and appropriate species.  Below LAYTROP, the water
1008! vapor self-continuum is interpolated (in temperature) separately. 
1009
1010      layreffr = 6
1011      laysolfr = laytrop
1012
1013! Lower atmosphere loop
1014      do lay = 1, laytrop
1015         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
1016            laysolfr = min(lay+1,laytrop)
1017         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(23) + 1
1018         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(23) + 1
1019         inds = indself(lay)
1020         indf = indfor(lay)
1021
1022         do ig = 1, ng23
1023            tauray = colmol(lay) * rayl(ig)
1024            taug(lay,ngs22+ig) = colh2o(lay) * &
1025                (givfac * (fac00(lay) * absa(ind0,ig) + &
1026                 fac10(lay) * absa(ind0+1,ig) + &
1027                 fac01(lay) * absa(ind1,ig) + &
1028                 fac11(lay) * absa(ind1+1,ig)) + &
1029                 selffac(lay) * (selfref(inds,ig) + &
1030                 selffrac(lay) * &
1031                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
1032                 forfac(lay) * (forref(indf,ig) + &
1033                 forfrac(lay) * &
1034                 (forref(indf+1,ig) - forref(indf,ig)))) 
1035!            ssa(lay,ngs22+ig) = tauray/taug(lay,ngs22+ig)
1036            if (lay .eq. laysolfr) sfluxzen(ngs22+ig) = sfluxref(ig) 
1037            taur(lay,ngs22+ig) = tauray
1038         enddo
1039      enddo
1040
1041! Upper atmosphere loop
1042      do lay = laytrop+1, nlayers
1043         do ig = 1, ng23
1044!            taug(lay,ngs22+ig) = colmol(lay) * rayl(ig)
1045!            ssa(lay,ngs22+ig) = 1.0_rb
1046            taug(lay,ngs22+ig) = 0._rb
1047            taur(lay,ngs22+ig) = colmol(lay) * rayl(ig) 
1048         enddo
1049      enddo
1050
1051      end subroutine taumol23
1052
1053!----------------------------------------------------------------------------
1054      subroutine taumol24
1055!----------------------------------------------------------------------------
1056!
1057!     band 24:  12850-16000 cm-1 (low - h2o,o2; high - o2)
1058!
1059!----------------------------------------------------------------------------
1060
1061! ------- Modules -------
1062
1063      use parrrsw, only : ng24, ngs23
1064      use rrsw_kg24, only : absa, ka, absb, kb, forref, selfref, &
1065                            sfluxref, abso3a, abso3b, rayla, raylb
1066
1067! ------- Declarations -------
1068
1069! Local
1070
1071      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
1072                          layreffr
1073      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1074                       fac110, fac111, fs, speccomb, specmult, specparm, &
1075                       tauray, strrat
1076
1077! Compute the optical depth by interpolating in ln(pressure),
1078! temperature, and appropriate species.  Below LAYTROP, the water
1079! vapor self-continuum is interpolated (in temperature) separately. 
1080
1081      strrat = 0.124692_rb
1082      layreffr = 1
1083      laysolfr = laytrop
1084
1085! Lower atmosphere loop
1086      do lay = 1, laytrop
1087         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
1088            laysolfr = min(lay+1,laytrop)
1089         speccomb = colh2o(lay) + strrat*colo2(lay)
1090         specparm = colh2o(lay)/speccomb 
1091         if (specparm .ge. oneminus) specparm = oneminus
1092         specmult = 8._rb*(specparm)
1093         js = 1 + int(specmult)
1094         fs = mod(specmult, 1._rb )
1095         fac000 = (1._rb - fs) * fac00(lay)
1096         fac010 = (1._rb - fs) * fac10(lay)
1097         fac100 = fs * fac00(lay)
1098         fac110 = fs * fac10(lay)
1099         fac001 = (1._rb - fs) * fac01(lay)
1100         fac011 = (1._rb - fs) * fac11(lay)
1101         fac101 = fs * fac01(lay)
1102         fac111 = fs * fac11(lay)
1103         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(24) + js
1104         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(24) + js
1105         inds = indself(lay)
1106         indf = indfor(lay)
1107
1108         do ig = 1, ng24
1109            tauray = colmol(lay) * (rayla(ig,js) + &
1110               fs * (rayla(ig,js+1) - rayla(ig,js)))
1111            taug(lay,ngs23+ig) = speccomb * &
1112                (fac000 * absa(ind0,ig) + &
1113                 fac100 * absa(ind0+1,ig) + &
1114                 fac010 * absa(ind0+9,ig) + &
1115                 fac110 * absa(ind0+10,ig) + &
1116                 fac001 * absa(ind1,ig) + &
1117                 fac101 * absa(ind1+1,ig) + &
1118                 fac011 * absa(ind1+9,ig) + &
1119                 fac111 * absa(ind1+10,ig)) + &
1120                 colo3(lay) * abso3a(ig) + &
1121                 colh2o(lay) * & 
1122                 (selffac(lay) * (selfref(inds,ig) + &
1123                 selffrac(lay) * &
1124                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
1125                 forfac(lay) * (forref(indf,ig) + & 
1126                 forfrac(lay) * &
1127                 (forref(indf+1,ig) - forref(indf,ig))))
1128!            ssa(lay,ngs23+ig) = tauray/taug(lay,ngs23+ig)
1129            if (lay .eq. laysolfr) sfluxzen(ngs23+ig) = sfluxref(ig,js) &
1130               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
1131            taur(lay,ngs23+ig) = tauray
1132         enddo
1133      enddo
1134
1135! Upper atmosphere loop
1136      do lay = laytrop+1, nlayers
1137         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(24) + 1
1138         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(24) + 1
1139
1140         do ig = 1, ng24
1141            tauray = colmol(lay) * raylb(ig)
1142            taug(lay,ngs23+ig) = colo2(lay) * &
1143                (fac00(lay) * absb(ind0,ig) + &
1144                 fac10(lay) * absb(ind0+1,ig) + &
1145                 fac01(lay) * absb(ind1,ig) + &
1146                 fac11(lay) * absb(ind1+1,ig)) + &
1147                 colo3(lay) * abso3b(ig)
1148!            ssa(lay,ngs23+ig) = tauray/taug(lay,ngs23+ig)
1149            taur(lay,ngs23+ig) = tauray
1150         enddo
1151      enddo
1152
1153      end subroutine taumol24
1154
1155!----------------------------------------------------------------------------
1156      subroutine taumol25
1157!----------------------------------------------------------------------------
1158!
1159!     band 25:  16000-22650 cm-1 (low - h2o; high - nothing)
1160!
1161!----------------------------------------------------------------------------
1162
1163! ------- Modules -------
1164
1165      use parrrsw, only : ng25, ngs24
1166      use rrsw_kg25, only : absa, ka, &
1167                            sfluxref, abso3a, abso3b, rayl
1168
1169! ------- Declarations -------
1170
1171! Local
1172
1173      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
1174                          layreffr
1175      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1176                       fac110, fac111, fs, speccomb, specmult, specparm, &
1177                       tauray
1178
1179! Compute the optical depth by interpolating in ln(pressure),
1180! temperature, and appropriate species.  Below LAYTROP, the water
1181! vapor self-continuum is interpolated (in temperature) separately. 
1182
1183      layreffr = 2
1184      laysolfr = laytrop
1185
1186! Lower atmosphere loop
1187      do lay = 1, laytrop
1188         if (jp(lay) .lt. layreffr .and. jp(lay+1) .ge. layreffr) &
1189            laysolfr = min(lay+1,laytrop)
1190         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(25) + 1
1191         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(25) + 1
1192
1193         do ig = 1, ng25
1194            tauray = colmol(lay) * rayl(ig)
1195            taug(lay,ngs24+ig) = colh2o(lay) * &
1196                (fac00(lay) * absa(ind0,ig) + &
1197                 fac10(lay) * absa(ind0+1,ig) + &
1198                 fac01(lay) * absa(ind1,ig) + &
1199                 fac11(lay) * absa(ind1+1,ig)) + &
1200                 colo3(lay) * abso3a(ig) 
1201!            ssa(lay,ngs24+ig) = tauray/taug(lay,ngs24+ig)
1202            if (lay .eq. laysolfr) sfluxzen(ngs24+ig) = sfluxref(ig) 
1203            taur(lay,ngs24+ig) = tauray
1204         enddo
1205      enddo
1206
1207! Upper atmosphere loop
1208      do lay = laytrop+1, nlayers
1209         do ig = 1, ng25
1210            tauray = colmol(lay) * rayl(ig)
1211            taug(lay,ngs24+ig) = colo3(lay) * abso3b(ig) 
1212!            ssa(lay,ngs24+ig) = tauray/taug(lay,ngs24+ig)
1213            taur(lay,ngs24+ig) = tauray
1214         enddo
1215      enddo
1216
1217      end subroutine taumol25
1218
1219!----------------------------------------------------------------------------
1220      subroutine taumol26
1221!----------------------------------------------------------------------------
1222!
1223!     band 26:  22650-29000 cm-1 (low - nothing; high - nothing)
1224!
1225!----------------------------------------------------------------------------
1226
1227! ------- Modules -------
1228
1229      use parrrsw, only : ng26, ngs25
1230      use rrsw_kg26, only : sfluxref, rayl
1231
1232! ------- Declarations -------
1233
1234! Local
1235
1236      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr
1237      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1238                       fac110, fac111, fs, speccomb, specmult, specparm, &
1239                       tauray
1240
1241! Compute the optical depth by interpolating in ln(pressure),
1242! temperature, and appropriate species.  Below LAYTROP, the water
1243! vapor self-continuum is interpolated (in temperature) separately. 
1244
1245      laysolfr = laytrop
1246
1247! Lower atmosphere loop
1248      do lay = 1, laytrop
1249         do ig = 1, ng26 
1250!            taug(lay,ngs25+ig) = colmol(lay) * rayl(ig)
1251!            ssa(lay,ngs25+ig) = 1.0_rb
1252            if (lay .eq. laysolfr) sfluxzen(ngs25+ig) = sfluxref(ig) 
1253            taug(lay,ngs25+ig) = 0._rb
1254            taur(lay,ngs25+ig) = colmol(lay) * rayl(ig) 
1255         enddo
1256      enddo
1257
1258! Upper atmosphere loop
1259      do lay = laytrop+1, nlayers
1260         do ig = 1, ng26
1261!            taug(lay,ngs25+ig) = colmol(lay) * rayl(ig)
1262!            ssa(lay,ngs25+ig) = 1.0_rb
1263            taug(lay,ngs25+ig) = 0._rb
1264            taur(lay,ngs25+ig) = colmol(lay) * rayl(ig) 
1265         enddo
1266      enddo
1267
1268      end subroutine taumol26
1269
1270!----------------------------------------------------------------------------
1271      subroutine taumol27
1272!----------------------------------------------------------------------------
1273!
1274!     band 27:  29000-38000 cm-1 (low - o3; high - o3)
1275!
1276!----------------------------------------------------------------------------
1277
1278! ------- Modules -------
1279
1280      use parrrsw, only : ng27, ngs26
1281      use rrsw_kg27, only : absa, ka, absb, kb, sfluxref, rayl
1282
1283! ------- Declarations -------
1284
1285! Local
1286
1287      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
1288                          layreffr
1289      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1290                       fac110, fac111, fs, speccomb, specmult, specparm, &
1291                       tauray, scalekur
1292
1293! Kurucz solar source function
1294! The values in sfluxref were obtained using the "low resolution"
1295! version of the Kurucz solar source function.  For unknown reasons,
1296! the total irradiance in this band differs from the corresponding
1297! total in the "high-resolution" version of the Kurucz function.
1298! Therefore, these values are scaled below by the factor SCALEKUR.
1299
1300      scalekur = 50.15_rb/48.37_rb
1301
1302! Compute the optical depth by interpolating in ln(pressure),
1303! temperature, and appropriate species.  Below LAYTROP, the water
1304! vapor self-continuum is interpolated (in temperature) separately. 
1305
1306      layreffr = 32
1307
1308! Lower atmosphere loop
1309      do lay = 1, laytrop
1310         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(27) + 1
1311         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(27) + 1
1312
1313         do ig = 1, ng27
1314            tauray = colmol(lay) * rayl(ig)
1315            taug(lay,ngs26+ig) = colo3(lay) * &
1316                (fac00(lay) * absa(ind0,ig) + &
1317                 fac10(lay) * absa(ind0+1,ig) + &
1318                 fac01(lay) * absa(ind1,ig) + &
1319                 fac11(lay) * absa(ind1+1,ig))
1320!            ssa(lay,ngs26+ig) = tauray/taug(lay,ngs26+ig)
1321            taur(lay,ngs26+ig) = tauray
1322         enddo
1323      enddo
1324
1325      laysolfr = nlayers
1326
1327! Upper atmosphere loop
1328      do lay = laytrop+1, nlayers
1329         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
1330            laysolfr = lay
1331         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(27) + 1
1332         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(27) + 1
1333
1334         do ig = 1, ng27
1335            tauray = colmol(lay) * rayl(ig)
1336            taug(lay,ngs26+ig) = colo3(lay) * &
1337                (fac00(lay) * absb(ind0,ig) + &
1338                 fac10(lay) * absb(ind0+1,ig) + &
1339                 fac01(lay) * absb(ind1,ig) + & 
1340                 fac11(lay) * absb(ind1+1,ig))
1341!            ssa(lay,ngs26+ig) = tauray/taug(lay,ngs26+ig)
1342            if (lay.eq.laysolfr) sfluxzen(ngs26+ig) = scalekur * sfluxref(ig) 
1343            taur(lay,ngs26+ig) = tauray
1344         enddo
1345      enddo
1346
1347      end subroutine taumol27
1348
1349!----------------------------------------------------------------------------
1350      subroutine taumol28
1351!----------------------------------------------------------------------------
1352!
1353!     band 28:  38000-50000 cm-1 (low - o3,o2; high - o3,o2)
1354!
1355!----------------------------------------------------------------------------
1356
1357! ------- Modules -------
1358
1359      use parrrsw, only : ng28, ngs27
1360      use rrsw_kg28, only : absa, ka, absb, kb, sfluxref, rayl
1361
1362! ------- Declarations -------
1363
1364! Local
1365
1366      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
1367                          layreffr
1368      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1369                       fac110, fac111, fs, speccomb, specmult, specparm, &
1370                       tauray, strrat
1371
1372! Compute the optical depth by interpolating in ln(pressure),
1373! temperature, and appropriate species.  Below LAYTROP, the water
1374! vapor self-continuum is interpolated (in temperature) separately. 
1375
1376      strrat = 6.67029e-07_rb
1377      layreffr = 58
1378
1379! Lower atmosphere loop
1380      do lay = 1, laytrop
1381         speccomb = colo3(lay) + strrat*colo2(lay)
1382         specparm = colo3(lay)/speccomb 
1383         if (specparm .ge. oneminus) specparm = oneminus
1384         specmult = 8._rb*(specparm)
1385         js = 1 + int(specmult)
1386         fs = mod(specmult, 1._rb )
1387         fac000 = (1._rb - fs) * fac00(lay)
1388         fac010 = (1._rb - fs) * fac10(lay)
1389         fac100 = fs * fac00(lay)
1390         fac110 = fs * fac10(lay)
1391         fac001 = (1._rb - fs) * fac01(lay)
1392         fac011 = (1._rb - fs) * fac11(lay)
1393         fac101 = fs * fac01(lay)
1394         fac111 = fs * fac11(lay)
1395         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(28) + js
1396         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(28) + js
1397         tauray = colmol(lay) * rayl
1398
1399         do ig = 1, ng28
1400            taug(lay,ngs27+ig) = speccomb * &
1401                (fac000 * absa(ind0,ig) + &
1402                 fac100 * absa(ind0+1,ig) + &
1403                 fac010 * absa(ind0+9,ig) + &
1404                 fac110 * absa(ind0+10,ig) + &
1405                 fac001 * absa(ind1,ig) + &
1406                 fac101 * absa(ind1+1,ig) + &
1407                 fac011 * absa(ind1+9,ig) + &
1408                 fac111 * absa(ind1+10,ig)) 
1409!            ssa(lay,ngs27+ig) = tauray/taug(lay,ngs27+ig)
1410            taur(lay,ngs27+ig) = tauray
1411         enddo
1412      enddo
1413
1414      laysolfr = nlayers
1415
1416! Upper atmosphere loop
1417      do lay = laytrop+1, nlayers
1418         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
1419            laysolfr = lay
1420         speccomb = colo3(lay) + strrat*colo2(lay)
1421         specparm = colo3(lay)/speccomb 
1422         if (specparm .ge. oneminus) specparm = oneminus
1423         specmult = 4._rb*(specparm)
1424         js = 1 + int(specmult)
1425         fs = mod(specmult, 1._rb )
1426         fac000 = (1._rb - fs) * fac00(lay)
1427         fac010 = (1._rb - fs) * fac10(lay)
1428         fac100 = fs * fac00(lay)
1429         fac110 = fs * fac10(lay)
1430         fac001 = (1._rb - fs) * fac01(lay)
1431         fac011 = (1._rb - fs) * fac11(lay)
1432         fac101 = fs * fac01(lay)
1433         fac111 = fs * fac11(lay)
1434         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(28) + js
1435         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(28) + js
1436         tauray = colmol(lay) * rayl
1437
1438         do ig = 1, ng28
1439            taug(lay,ngs27+ig) = speccomb * &
1440                (fac000 * absb(ind0,ig) + &
1441                 fac100 * absb(ind0+1,ig) + &
1442                 fac010 * absb(ind0+5,ig) + &
1443                 fac110 * absb(ind0+6,ig) + &
1444                 fac001 * absb(ind1,ig) + &
1445                 fac101 * absb(ind1+1,ig) + &
1446                 fac011 * absb(ind1+5,ig) + &
1447                 fac111 * absb(ind1+6,ig)) 
1448!            ssa(lay,ngs27+ig) = tauray/taug(lay,ngs27+ig)
1449            if (lay .eq. laysolfr) sfluxzen(ngs27+ig) = sfluxref(ig,js) &
1450               + fs * (sfluxref(ig,js+1) - sfluxref(ig,js))
1451            taur(lay,ngs27+ig) = tauray
1452         enddo
1453      enddo
1454
1455      end subroutine taumol28
1456
1457!----------------------------------------------------------------------------
1458      subroutine taumol29
1459!----------------------------------------------------------------------------
1460!
1461!     band 29:  820-2600 cm-1 (low - h2o; high - co2)
1462!
1463!----------------------------------------------------------------------------
1464
1465! ------- Modules -------
1466
1467      use parrrsw, only : ng29, ngs28
1468      use rrsw_kg29, only : absa, ka, absb, kb, forref, selfref, &
1469                            sfluxref, absh2o, absco2, rayl
1470
1471! ------- Declarations -------
1472
1473! Local
1474
1475      integer(kind=im) :: ig, ind0, ind1, inds, indf, js, lay, laysolfr, &
1476                          layreffr
1477      real(kind=rb) :: fac000, fac001, fac010, fac011, fac100, fac101, &
1478                       fac110, fac111, fs, speccomb, specmult, specparm, &
1479                       tauray
1480
1481! Compute the optical depth by interpolating in ln(pressure),
1482! temperature, and appropriate species.  Below LAYTROP, the water
1483! vapor self-continuum is interpolated (in temperature) separately. 
1484
1485      layreffr = 49
1486
1487! Lower atmosphere loop
1488      do lay = 1, laytrop
1489         ind0 = ((jp(lay)-1)*5+(jt(lay)-1))*nspa(29) + 1
1490         ind1 = (jp(lay)*5+(jt1(lay)-1))*nspa(29) + 1
1491         inds = indself(lay)
1492         indf = indfor(lay)
1493         tauray = colmol(lay) * rayl
1494
1495         do ig = 1, ng29
1496            taug(lay,ngs28+ig) = colh2o(lay) * &
1497               ((fac00(lay) * absa(ind0,ig) + &
1498                 fac10(lay) * absa(ind0+1,ig) + &
1499                 fac01(lay) * absa(ind1,ig) + &
1500                 fac11(lay) * absa(ind1+1,ig)) + &
1501                 selffac(lay) * (selfref(inds,ig) + &
1502                 selffrac(lay) * &
1503                 (selfref(inds+1,ig) - selfref(inds,ig))) + &
1504                 forfac(lay) * (forref(indf,ig) + & 
1505                 forfrac(lay) * &
1506                 (forref(indf+1,ig) - forref(indf,ig)))) &
1507                 + colco2(lay) * absco2(ig) 
1508!            ssa(lay,ngs28+ig) = tauray/taug(lay,ngs28+ig)
1509            taur(lay,ngs28+ig) = tauray
1510         enddo
1511      enddo
1512
1513      laysolfr = nlayers
1514
1515! Upper atmosphere loop
1516      do lay = laytrop+1, nlayers
1517         if (jp(lay-1) .lt. layreffr .and. jp(lay) .ge. layreffr) &
1518            laysolfr = lay
1519         ind0 = ((jp(lay)-13)*5+(jt(lay)-1))*nspb(29) + 1
1520         ind1 = ((jp(lay)-12)*5+(jt1(lay)-1))*nspb(29) + 1
1521         tauray = colmol(lay) * rayl
1522
1523         do ig = 1, ng29
1524            taug(lay,ngs28+ig) = colco2(lay) * &
1525                (fac00(lay) * absb(ind0,ig) + &
1526                 fac10(lay) * absb(ind0+1,ig) + &
1527                 fac01(lay) * absb(ind1,ig) + &
1528                 fac11(lay) * absb(ind1+1,ig)) & 
1529                 + colh2o(lay) * absh2o(ig) 
1530!            ssa(lay,ngs28+ig) = tauray/taug(lay,ngs28+ig)
1531            if (lay .eq. laysolfr) sfluxzen(ngs28+ig) = sfluxref(ig) 
1532            taur(lay,ngs28+ig) = tauray
1533         enddo
1534      enddo
1535
1536      end subroutine taumol29
1537
1538      end subroutine taumol_sw
1539
1540      end module rrtmg_sw_taumol
1541
Note: See TracBrowser for help on using the repository browser.