emm(i,2)    = 0.5*(co2em(i,k2) + co2eml(i,k2))               tbar(i,3)   = 0.5*(tbar(i,2) + tbar(i,1))               emm(i,3)    = emm(i,1)               tbar(i,4)   = tbar(i,3)               emm(i,4)    = emm(i,2)               o3emm(i,1)  = 0.5*(dbvtit(i,k2+1) + dbvtly(i,k2))               o3emm(i,2)  = 0.5*(dbvtit(i,k2) + dbvtly(i,k2))               o3emm(i,3)  = o3emm(i,1)               o3emm(i,4)  = o3emm(i,2)               temh2o(i,1) = tbar(i,1)               temh2o(i,2) = tbar(i,2)               temh2o(i,3) = tbar(i,1)               temh2o(i,4) = tbar(i,2)               dpnm(i)     = pnm(i,k2+1) - pnm(i,k2)            end do!!  Weighted Planck functions for trace gases!            do wvl = 1,14               do i = 1,ncol                  bplnk(wvl,i,1) = 0.5*(abplnk1(wvl,i,k2+1) + abplnk2(wvl,i,k2))                  bplnk(wvl,i,2) = 0.5*(abplnk1(wvl,i,k2) + abplnk2(wvl,i,k2))                  bplnk(wvl,i,3) = bplnk(wvl,i,1)                  bplnk(wvl,i,4) = bplnk(wvl,i,2)               end do            end do            do i=1,ncol               rdpnmsq    = 1./(pnmsq(i,k2+1) - pnmsq(i,k2))               rdpnm      = 1./dpnm(i)               p1         = .5*(pbr(i,k2) + pnm(i,k2+1))               p2         = .5*(pbr(i,k2) + pnm(i,k2  ))               uinpl(i,1) =  (pnmsq(i,k2+1) - p1**2)*rdpnmsq               uinpl(i,2) = -(pnmsq(i,k2  ) - p2**2)*rdpnmsq               uinpl(i,3) = -(pnmsq(i,k2  ) - p1**2)*rdpnmsq               uinpl(i,4) =  (pnmsq(i,k2+1) - p2**2)*rdpnmsq               winpl(i,1) = (.5*( pnm(i,k2+1) - pbr(i,k2)))*rdpnm               winpl(i,2) = (.5*(-pnm(i,k2  ) + pbr(i,k2)))*rdpnm               winpl(i,3) = (.5*( pnm(i,k2+1) + pbr(i,k2)) - pnm(i,k2  ))*rdpnm               winpl(i,4) = (.5*(-pnm(i,k2  ) - pbr(i,k2)) + pnm(i,k2+1))*rdpnm               tmp1       = 1./(piln(i,k2+1) - piln(i,k2))               tmp2       = piln(i,k2+1) - pmln(i,k2)               tmp3       = piln(i,k2  ) - pmln(i,k2)               zinpl(i,1) = (.5*tmp2          )*tmp1               zinpl(i,2) = (        - .5*tmp3)*tmp1               zinpl(i,3) = (.5*tmp2 -    tmp3)*tmp1               zinpl(i,4) = (   tmp2 - .5*tmp3)*tmp1               pinpl(i,1) = 0.5*(p1 + pnm(i,k2+1))               pinpl(i,2) = 0.5*(p2 + pnm(i,k2  ))               pinpl(i,3) = 0.5*(p1 + pnm(i,k2  ))               pinpl(i,4) = 0.5*(p2 + pnm(i,k2+1))            end do            do 400 kn=1,4               do i=1,ncol                  u(i)     = uinpl(i,kn)*abs(plh2o(i,k2) - plh2o(i,k2+1))                  sqrtu(i) = sqrt(u(i))                  dw(i)    = abs(w(i,k2) - w(i,k2+1))                  pnew(i)  = u(i)/(winpl(i,kn)*dw(i))                  pnew_mks  = pnew(i) * sslp_mks                  t_p = min(max(tbar(i,kn), 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)                  q_path = dw(i) / ABS(dpnm(i)) / rga                  ds2c     = abs(s2c(i,k2) - s2c(i,k2+1))                  uc1(i)   = uinpl(i,kn)*ds2c                  pch2o    = uc1(i)                  uc1(i)   = (uc1(i) + 1.7e-3*u(i))*(1. +  2.*uc1(i))/(1. + 15.*uc1(i))                  dtx(i)      = temh2o(i,kn) - 250.                  dty(i)      = tbar(i,kn) - 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  = temh2o(i,kn)                  te2  = te1 * te1                  te3  = te2 * te1                  te4  = te3 * te1                  te5  = te4 * te1!! Indices for lines and continuum tables ! Note: because we are dealing with the nearest layer,!       the Hulst-Curtis-Godson corrections!       for inhomogeneous paths are not applied.!                  uvar = u(i)*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(pnew(i), 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                                    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(temh2o(i,kn)-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                                          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 !! Non-window absorptivity!                  ib = 1                  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!                  if (uvar < min_u_h2o) then                     uscl = uvar / min_u_h2o                     abso(i,ib) = abso(i,ib) * uscl                  endif                               !! Window absorptivity!                  ib = 2                  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 = &                       ah2ow(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &                       ah2ow(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &                       ah2ow(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &                       ah2ow(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &                       ah2ow(ip , itp , iu1, ite , irh ) * w11_11 * wu  + &                       ah2ow(ip , itp , iu1, ite , irh1) * w11_10 * wu  + &                       ah2ow(ip , itp , iu1, ite1, irh ) * w11_01 * wu  + &                       ah2ow(ip , itp , iu1, ite1, irh1) * w11_00 * wu  + &                       ah2ow(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &                       ah2ow(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &                       ah2ow(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &                       ah2ow(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &                       ah2ow(ip , itp1, iu1, ite , irh ) * w10_11 * wu  + &                       ah2ow(ip , itp1, iu1, ite , irh1) * w10_10 * wu  + &                       ah2ow(ip , itp1, iu1, ite1, irh ) * w10_01 * wu  + &                       ah2ow(ip , itp1, iu1, ite1, irh1) * w10_00 * wu  + &                       ah2ow(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &                       ah2ow(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &                       ah2ow(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &                       ah2ow(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &                       ah2ow(ip1, itp , iu1, ite , irh ) * w01_11 * wu  + &                       ah2ow(ip1, itp , iu1, ite , irh1) * w01_10 * wu  + &                       ah2ow(ip1, itp , iu1, ite1, irh ) * w01_01 * wu  + &                       ah2ow(ip1, itp , iu1, ite1, irh1) * w01_00 * wu  + &                       ah2ow(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &                       ah2ow(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &                       ah2ow(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &                       ah2ow(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &                       ah2ow(ip1, itp1, iu1, ite , irh ) * w00_11 * wu  + &                       ah2ow(ip1, itp1, iu1, ite , irh1) * w00_10 * wu  + &                       ah2ow(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu  + &                       ah2ow(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!                  if (uvar < min_u_