Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

rrtmg_sw_taumol.f90 @ 4559

Last change on this file since 4559 was 1585, checked in by maronga, 10 years ago
Added support for RRTMG radiation code
File size: 61.9 KB

Rev	Line
[1585]	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(65NSPA(n),MG), ABSB(235NSPB(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_rbcolo2(lay)/(350.0_rb2.0_rb)
	912	speccomb = colh2o(lay) + o2adjstrratcolo2(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_rbcolo2(lay)/(350.0_rb2.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.

Context Navigation

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

Download in other formats: