r250    = 1./250.   r3205   = 1./.3205   r300    = 1./300.   rsslp   = 1./sslp   r2sslp  = 1./(2.*sslp)!!Constants for computing U corresponding to H2O cont. path!   fdif       = 1.66   sslp_mks   = sslp / 10.0   omeps      = 1.0 - epsilo!! Non-adjacent layer absorptivity:!! abso(i,1)     0 -  800 cm-1   h2o rotation band! abso(i,1)  1200 - 2200 cm-1   h2o vibration-rotation band! abso(i,2)   800 - 1200 cm-1   h2o window!! Separation between rotation and vibration-rotation dropped, so!                only 2 slots needed for H2O absorptivity!! 500-800 cm^-1 H2o continuum/line overlap already included!                in abso(i,1).  This used to be in abso(i,4)!! abso(i,3)   o3  9.6 micrometer band (nu3 and nu1 bands)! abso(i,4)   co2 15  micrometer band system!   do k=ntoplw,pverp      do i=1,ncol         pnmsq(i,k) = pnm(i,k)**2         dtx(i) = tplnka(i,k) - 250.      end do   end do!! Non-nearest layer level loops!   do 200 k1=pverp,ntoplw,-1      do 100 k2=pverp,ntoplw,-1         if (k1 == k2) go to 100         do i=1,ncol            dplh2o(i) = plh2o(i,k1) - plh2o(i,k2)            u(i)      = abs(dplh2o(i))            sqrtu(i)  = sqrt(u(i))            ds2c      = abs(s2c(i,k1) - s2c(i,k2))            dw(i)     = abs(w(i,k1) - w(i,k2))            uc1(i)    = (ds2c + 1.7e-3*u(i))*(1. +  2.*ds2c)/(1. + 15.*ds2c)            pch2o     = ds2c            pnew(i)   = u(i)/dw(i)            pnew_mks  = pnew(i) * sslp_mks!! Changed effective path temperature to std. Curtis-Godson form!            tpatha = abs(tcg(i,k1) - tcg(i,k2))/dw(i)            t_p = min(max(tpatha, min_tp_h2o), max_tp_h2o)            iest = floor(t_p) - min_tp_h2o            esx = estblh2o(iest) + (estblh2o(iest+1)-estblh2o(iest)) * &                 (t_p - min_tp_h2o - iest)            qsx = epsilo * esx / (pnew_mks - omeps * esx)!! Compute effective RH along path!            q_path = dw(i) / abs(pnm(i,k1) - pnm(i,k2)) / rga!! Calculate effective u, pnew for each band using!        Hulst-Curtis-Godson approximation:! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis, !           2nd edition, Oxford University Press, 1989.! Effective H2O path (w)!      eq. 6.24, p. 228! Effective H2O path pressure (pnew = u/w):!      eq. 6.29, p. 228!            ub(1) = abs(plh2ob(1,i,k1) - plh2ob(1,i,k2)) / psi(t_p,1)            ub(2) = abs(plh2ob(2,i,k1) - plh2ob(2,i,k2)) / psi(t_p,2)                        pnewb(1) = ub(1) / abs(wb(1,i,k1) - wb(1,i,k2)) * phi(t_p,1)            pnewb(2) = ub(2) / abs(wb(2,i,k1) - wb(2,i,k2)) * phi(t_p,2)            dtx(i)      = tplnka(i,k2) - 250.            dty(i)      = tpatha       - 250.            fwk(i)  = fwcoef + fwc1/(1. + fwc2*u(i))            fwku(i) = fwk(i)*u(i)!! Define variables for C/H/E fit!! abso(i,1)     0 -  800 cm-1   h2o rotation band! abso(i,1)  1200 - 2200 cm-1   h2o vibration-rotation band! abso(i,2)   800 - 1200 cm-1   h2o window!! Separation between rotation and vibration-rotation dropped, so!                only 2 slots needed for H2O absorptivity!! Notation:! U   = integral (P/P_0 dW)  ! P   = atmospheric pressure! P_0 = reference atmospheric pressure! W   = precipitable water path! T_e = emission temperature! T_p = path temperature! RH  = path relative humidity!!! Terms for asymptotic value of emissivity!            te1  = tplnka(i,k2)            te2  = te1 * te1            te3  = te2 * te1            te4  = te3 * te1            te5  = te4 * te1!!  Band-independent indices for lines and continuum tables!            dvar = (t_p - min_tp_h2o) / dtp_h2o            itp = min(max(int(aint(dvar,r8)) + 1, 1), n_tp - 1)            itp1 = itp + 1            wtp = dvar - floor(dvar)            wtp1 = 1.0 - wtp                        t_e = min(max(tplnka(i,k2)-t_p, min_te_h2o), max_te_h2o)            dvar = (t_e - min_te_h2o) / dte_h2o            ite = min(max(int(aint(dvar,r8)) + 1, 1), n_te - 1)            ite1 = ite + 1            wte = dvar - floor(dvar)            wte1 = 1.0 - wte                        rh_path = min(max(q_path / qsx, min_rh_h2o), max_rh_h2o)            dvar = (rh_path - min_rh_h2o) / drh_h2o            irh = min(max(int(aint(dvar,r8)) + 1, 1), n_rh - 1)            irh1 = irh + 1            wrh = dvar - floor(dvar)            wrh1 = 1.0 - wrh            w_0_0_ = wtp  * wte            w_0_1_ = wtp  * wte1            w_1_0_ = wtp1 * wte             w_1_1_ = wtp1 * wte1                        w_0_00 = w_0_0_ * wrh            w_0_01 = w_0_0_ * wrh1            w_0_10 = w_0_1_ * wrh            w_0_11 = w_0_1_ * wrh1            w_1_00 = w_1_0_ * wrh            w_1_01 = w_1_0_ * wrh1            w_1_10 = w_1_1_ * wrh            w_1_11 = w_1_1_ * wrh1!! H2O Continuum path for 0-800 and 1200-2200 cm^-1!!    Assume foreign continuum dominates total H2O continuum in these bands!    per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776!    Then the effective H2O path is just !         U_c = integral[ f(P) dW ]!    where !           W = water-vapor mass and !        f(P) = dependence of foreign continuum on pressure !             = P / sslp!    Then !         U_c = U (the same effective H2O path as for lines)!!! Continuum terms for 800-1200 cm^-1!!    Assume self continuum dominates total H2O continuum for this band!    per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776!    Then the effective H2O self-continuum path is !         U_c = integral[ h(e,T) dW ]                        (*eq. 1*)!    where !           W = water-vapor mass and !           e = partial pressure of H2O along path!           T = temperature along path!      h(e,T) = dependence of foreign continuum on e,T!             = e / sslp * f(T)!!    Replacing!           e =~ q * P / epsilo!           q = mixing ratio of H2O!     epsilo = 0.622!!    and using the definition!           U = integral [ (P / sslp) dW ]!             = (P / sslp) W                                 (homogeneous path)!!    the effective path length for the self continuum is!         U_c = (q / epsilo) f(T) U                         (*eq. 2*)!!    Once values of T, U, and q have been calculated for the inhomogeneous!        path, this sets U_c for the corresponding!        homogeneous atmosphere.  However, this need not equal the!        value of U_c' defined by eq. 1 for the actual inhomogeneous atmosphere!        under consideration.!!    Solution: hold T and q constant, solve for U' that gives U_c' by!        inverting eq. (2):!!        U' = (U_c * epsilo) / (q * f(T))!            fch2o = fh2oself(t_p)             uch2o = (pch2o * epsilo) / (q_path * fch2o)!! Band-dependent indices for non-window!            ib = 1            uvar = ub(ib) * fdif            log_u  = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)            dvar = (log_u - min_lu_h2o) / dlu_h2o            iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)            iu1 = iu + 1            wu = dvar - floor(dvar)            wu1 = 1.0 - wu                        log_p  = min(log10(max(pnewb(ib), min_p_h2o)), max_lp_h2o)            dvar = (log_p - min_lp_h2o) / dlp_h2o            ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)            ip1 = ip + 1            wp = dvar - floor(dvar)            wp1 = 1.0 - wp                     w00_00 = wp  * w_0_00             w00_01 = wp  * w_0_01             w00_10 = wp  * w_0_10             w00_11 = wp  * w_0_11             w01_00 = wp  * w_1_00             w01_01 = wp  * w_1_01             w01_10 = wp  * w_1_10             w01_11 = wp  * w_1_11             w10_00 = wp1 * w_0_00             w10_01 = wp1 * w_0_01             w10_10 = wp1 * w_0_10             w10_11 = wp1 * w_0_11             w11_00 = wp1 * w_1_00             w11_01 = wp1 * w_1_01             w11_10 = wp1 * w_1_10             w11_11 = wp1 * w_1_11 !! Asymptotic value of absorptivity as U->infinity!            fa = fat(1,ib) + &                 fat(2,ib) * te1 + &                 fat(3,ib) * te2 + &                 fat(4,ib) * te3 + &                 fat(5,ib) * te4 + &                 fat(6,ib) * te5            a_star = &                 ah2onw(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &                 ah2onw(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &                 ah2onw(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &                 ah2onw(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &                 ah2onw(ip , itp , iu1, ite , irh ) * w11_11 * wu  + &                 ah2onw(ip , itp , iu1, ite , irh1) * w11_10 * wu  + &                 ah2onw(ip , itp , iu1, ite1, irh ) * w11_01 * wu  + &                 ah2onw(ip , itp , iu1, ite1, irh1) * w11_00 * wu  + &                 ah2onw(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &                 ah2onw(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &                 ah2onw(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &                 ah2onw(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &                 ah2onw(ip , itp1, iu1, ite , irh ) * w10_11 * wu  + &                 ah2onw(ip , itp1, iu1, ite , irh1) * w10_10 * wu  + &                 ah2onw(ip , itp1, iu1, ite1, irh ) * w10_01 * wu  + &                 ah2onw(ip , itp1, iu1, ite1, irh1) * w10_00 * wu  + &                 ah2onw(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &                 ah2onw(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &                 ah2onw(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &                 ah2onw(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &                 ah2onw(ip1, itp , iu1, ite , irh ) * w01_11 * wu  + &                 ah2onw(ip1, itp , iu1, ite , irh1) * w01_10 * wu  + &                 ah2onw(ip1, itp , iu1, ite1, irh ) * w01_01 * wu  + &                 ah2onw(ip1, itp , iu1, ite1, irh1) * w01_00 * wu  + &                 ah2onw(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &                 ah2onw(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &                 ah2onw(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &                 ah2onw(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &                 ah2onw(ip1, itp1, iu1, ite , irh ) * w00_11 * wu  + &                 ah2onw(ip1, itp1, iu1, ite , irh1) * w00_10 * wu  + &                 ah2onw(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu  + &                 ah2onw(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu             abso(i,ib) = min(max(fa * a_star, 0.0_r8), 1.0_r8)!! Invoke linear limit for scaling wrt u below min_u_h2o