var/da/da_ssmi/da_tbatmos_tl.inc

   1 subroutine da_tbatmos_tl(ifreq,theta,p0,wv,hwv,ta,gamma,lw,zcld,   &
   2    tbup,tbdn,tauatm,                         &
   3    TGL_theta,TGL_p0,TGL_wv,TGL_hwv,TGL_ta,TGL_gamma,   &
   4    TGL_lw,TGL_zcld,TGL_tbup,TGL_tbdn,        &
   5    TGL_tauatm                               )
   6
   7    !-----------------------------------------------------------------
   8    ! Purpose: TBD
   9    ! Input  : TGL_p0,TGL_wv,TGL_hwv,TGL_ta,TGL_gamma,TGL_lw,TGL_zcld
  10    !          TGL_theta (somtime theta is a variable)
  11    ! Output : TGL_tbup,TGL_tbdn,TGL_tauatm,tbup,tbdn,tauatm
  12    ! -----------------------------------------------------------------
  13
  14    implicit none
  15
  16    integer,intent(in)  :: ifreq
  17    real   ,intent(in)  :: theta,p0,wv,hwv,ta,gamma,lw,zcld
  18    real   ,intent(in)  :: TGL_p0,TGL_wv,TGL_hwv,TGL_ta,TGL_gamma,TGL_lw,TGL_zcld,TGL_theta
  19    real   ,intent(out) :: tbup,tbdn,tauatm,TGL_tbup,TGL_tbdn, TGL_tauatm
  20
  21    real :: mu,hdn,hup,hdninf,hupinf,TGL_mu
  22
  23    real :: b1(4),b2(4),b3(4)
  24    real :: c(4),d1(4),d2(4),d3(4),zeta(4),kw0(4),kw1(4),kw2(4),kw3(4)
  25    real :: tau,tau1,tau2,taucld
  26    real :: tcld,tc,em,em1
  27    real :: sigv,sigo,sig,sig1,sigcld
  28    real :: teff1dn,teff1up,teffdn,teffup
  29    real :: tbcld,tbclrdn,tbclrup,tb1dn,tb1up,tb2dn,tb2up
  30    real :: otbar,tc2,tc3,hv,ho,alph
  31    real :: TGL_sigv,TGL_otbar,TGL_sigo,TGL_tcld,TGL_tc,TGL_tc2,TGL_tc3
  32    real :: TGL_sigcld,TGL_taucld,TGL_tbcld,TGL_hv,TGL_ho
  33    real :: TGL_hdn,TGL_hup,TGL_hdninf,TGL_sig,TGL_sig1,TGL_tau,TGL_tau1
  34    real :: TGL_tau2,TGL_em1,TGL_teff1dn,TGL_hupinf,TGL_em,TGL_teff1up
  35    real :: TGL_teffdn,TGL_teffup,TGL_tbclrdn,TGL_tbclrup,TGL_tb1dn,TGL_tb1up
  36    real :: TGL_tb2dn,TGL_tb2up,TGL_alph
  37
  38    data b1/-.46847e-1,-.57752e-1,-.18885,.10990/
  39    data b2/.26640e-4,.31662e-4,.9832e-4,.60531e-4/
  40    data b3/.87560e+1,.10961e+2,.36678e+2,-.37578e+2/
  41    data c/ .9207,   1.208,     .8253,     .8203/
  42    data zeta/4.2,4.2,4.2,2.9/
  43    data d1/-.35908e+1,-.38921e+1,-.43072e+1,-.17020e+0/
  44    data d2/ .29797e-1, .31054e-1, .32801e-1, .13610e-1/
  45    data d3/-.23174e-1,-.23543e-1,-.24101e-1,-.15776e+0/
  46    data kw0/ .786e-1, .103,    .267,    .988/
  47    data kw1/-.230e-2,-.296e-2,-.673e-2,-.107e-1/
  48    data kw2/ .448e-4, .557e-4, .975e-4,-.535e-4/
  49    data kw3/-.464e-6,-.558e-6,-.724e-6, .115e-5/
  50
  51    if (trace_use) call da_trace_entry("da_tbatmos_tl")
  52
  53    ! mu = secant(theta)
  54    ! somtime theta is a variable
  55
  56    mu     = 1.0/cos(theta*0.0174533)
  57    TGL_mu = mu*mu*0.0174533*TGL_theta*sin(theta*0.0174533)
  58
  59    ! get water vapor optical depth
  60
  61    call da_sigma_v_tl(ifreq,p0,wv,hwv,ta,gamma,sigv,   &
  62                         TGL_p0,TGL_wv,TGL_hwv,TGL_ta,    &
  63                         TGL_gamma,TGL_sigv              )
  64
  65    ! otbar = one over "mean" temperature
  66
  67        otbar =   1.0/(ta - gamma*zeta(ifreq))
  68    TGL_otbar = - otbar*otbar*(TGL_ta - TGL_gamma*zeta(ifreq))
  69
  70    ! sigo = dry air optical depth
  71
  72        sigo = b1(ifreq) + b2(ifreq)*    p0  + b3(ifreq)*    otbar
  73    TGL_sigo =             b2(ifreq)*TGL_p0  + b3(ifreq)*TGL_otbar
  74
  75    ! cloud parameters
  76
  77        tcld   =     ta - gamma*zcld
  78    TGL_tcld   = TGL_ta - TGL_gamma*zcld - gamma*TGL_zcld
  79          tc   =     tcld - t_kelvin
  80      TGL_tc   = TGL_tcld
  81         tc2   = tc*tc
  82     TGL_tc2   = 2.0*tc*TGL_tc
  83         tc3   = tc2*tc
  84     TGL_tc3   = TGL_tc2*tc + tc2*TGL_tc
  85        sigcld =  (kw0(ifreq) + tc*kw1(ifreq) + tc2*kw2(ifreq) +  &
  86                    tc3*kw3(ifreq))*lw
  87    TGL_sigcld =  (TGL_tc *kw1(ifreq) + TGL_tc2*kw2(ifreq) +  &
  88                    TGL_tc3*kw3(ifreq))*lw &
  89                + (kw0(ifreq) + tc*kw1(ifreq) + tc2*kw2(ifreq) +  &
  90                    tc3*kw3(ifreq))*TGL_lw
  91
  92        taucld =   exp(-mu*sigcld)
  93    TGL_taucld = - (TGL_mu*sigcld + mu*TGL_sigcld)*taucld
  94         tbcld =   (1.0 - taucld)*tcld
  95     TGL_tbcld = - TGL_taucld*tcld + (1.0 - taucld)*TGL_tcld
  96
  97    ! hv, ho = effective absorber scale heights for vapor, dry air
  98
  99        hv = c(ifreq)*    hwv
 100    TGL_hv = c(ifreq)*TGL_hwv
 101        ho = d1(ifreq) + d2(ifreq)*    ta + d3(ifreq)*    gamma
 102    TGL_ho =             d2(ifreq)*TGL_ta + d3(ifreq)*TGL_gamma
 103
 104    ! get effective emission heights for layer 1 and total atmosphere
 105
 106    call da_effht_tl(ho,hv,sigo,sigv,mu,zcld,hdn,hup,        &
 107                   hdninf,hupinf,                          &
 108                   TGL_ho,TGL_hv,TGL_sigo,TGL_sigv,TGL_mu, &
 109                   TGL_zcld,TGL_hdn,TGL_hup,TGL_hdninf,    &
 110                   TGL_hupinf                             )
 111
 112    ! atmospheric transmittances in layer one and two, and combined
 113
 114         sig =     sigo +     sigv
 115
 116     TGL_sig = TGL_sigo + TGL_sigv
 117        sig1 = sigo*(1.0-exp(-zcld/ho)) + sigv*(1.0-exp(-zcld/hv))
 118    TGL_sig1 =   TGL_sigo*(1.0-exp(-zcld/ho))  &
 119               + TGL_sigv*(1.0-exp(-zcld/hv))  &
 120               + sigo*(TGL_zcld/ho - zcld*TGL_ho/(ho*ho))*exp(-zcld/ho) &
 121               + sigv*(TGL_zcld/hv - zcld*TGL_hv/(hv*hv))*exp(-zcld/hv)
 122        tau  =  exp(-mu*sig)
 123    TGL_tau  = -(TGL_mu*sig + mu*TGL_sig) * tau
 124        tau1 =  exp(-mu*sig1)
 125    TGL_tau1 = -(mu*TGL_sig1 + TGL_mu*sig1) *tau1
 126        tau2 =  tau/tau1
 127    TGL_tau2 =  TGL_tau/tau1 - tau2*TGL_tau1/tau1
 128
 129    ! atmospheric "emissivity"
 130
 131        em1  =   1.0 - tau1
 132    TGL_em1  = - TGL_tau1
 133        em   =   1.0 - tau
 134    TGL_em   = - TGL_tau
 135
 136    ! downwelling and upwelling brightness temperature for each layer
 137
 138        teff1dn =     ta - gamma*hdn
 139    TGL_teff1dn = TGL_ta - TGL_gamma*hdn - gamma*TGL_hdn
 140        teff1up =     ta - gamma*hup
 141    TGL_teff1up = TGL_ta - TGL_gamma*hup - gamma*TGL_hup
 142         teffdn =     ta - gamma*hdninf
 143     TGL_teffdn = TGL_ta - TGL_gamma*hdninf - gamma*TGL_hdninf
 144         teffup =     ta - gamma*hupinf
 145     TGL_teffup = TGL_ta - TGL_gamma*hupinf - gamma*TGL_hupinf
 146        tbclrdn = teffdn*em
 147    TGL_tbclrdn = TGL_teffdn*em + teffdn*TGL_em
 148        tbclrup = teffup*em
 149    TGL_tbclrup = TGL_teffup*em + teffup*TGL_em
 150
 151         tb1dn = em1*teff1dn
 152    TGL_tb1dn = TGL_em1*teff1dn + em1*TGL_teff1dn
 153         tb1up = em1*teff1up
 154    TGL_tb1up = TGL_em1*teff1up + em1*TGL_teff1up
 155         tb2dn = (tbclrdn - tb1dn)/tau1
 156    TGL_tb2dn = (TGL_tbclrdn - TGL_tb1dn)/tau1 - tb2dn*TGL_tau1/tau1
 157         tb2up =      tbclrup - tau2*tb1up
 158    TGL_tb2up =  TGL_tbclrup - TGL_tau2*tb1up - tau2*TGL_tb1up
 159
 160    ! total downwelling and upwelling brightness temperature and transmittance
 161
 162        tbdn  =     tb1dn + tau1*(tbcld + taucld*tb2dn)
 163    TGL_tbdn  = TGL_tb1dn + TGL_tau1*(tbcld + taucld*tb2dn) &
 164                          + tau1*(TGL_tbcld + TGL_taucld*tb2dn + taucld*TGL_tb2dn)
 165        tbup  =     tb2up + tau2*(tbcld + taucld*tb1up)
 166    TGL_tbup  = TGL_tb2up + TGL_tau2*(tbcld + taucld*tb1up) &
 167                          + tau2*(TGL_tbcld + TGL_taucld*tb1up + taucld*TGL_tb1up)
 168        tauatm = tau*taucld
 169    TGL_tauatm = TGL_tau*taucld + tau*TGL_taucld
 170    !
 171    ! the following lines apply an ad hoc correction to improve fit
 172    ! at large angles and/or high gaseous opacities
 173    !  (downwelling brightness temperatures only)
 174
 175    alph = (0.636619*atan(mu*sig))**2
 176    if (abs(sig) .gt. 0.0) then
 177       TGL_alph = 2.0*0.636619*0.636619* &
 178                  (TGL_mu*sig + mu*TGL_sig)*atan(mu*sig)/(1.0+mu*mu*sig*sig)
 179    else
 180       TGL_alph = 0.0
 181    end if
 182    TGL_tbdn = - TGL_alph*tbdn + (1.0-alph)*TGL_tbdn &
 183               + TGL_em*alph*ta + em*TGL_alph*ta + em*alph*TGL_ta
 184    tbdn = (1.0-alph)*tbdn + em*alph*ta
 185
 186    if (trace_use) call da_trace_exit("da_tbatmos_tl")
 187
 188 end subroutine da_tbatmos_tl
 189
 190