var/da/da_ssmi/da_tbatmos_adj.inc

   1 subroutine da_tbatmos_adj(ifreq,theta,p0,wv,hwv,ta,gamma,lw,zcld,   &
   2    tbup,tbdn,tauatm, ADJ_theta,ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta,ADJ_gamma,   &
   3    ADJ_lw,ADJ_zcld,ADJ_tbup,ADJ_tbdn, ADJ_tauatm)
   4
   5    implicit none
   6
   7    !-----------------------------------------------------------------
   8    ! Purpose: TBD
   9    ! Output : ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta,ADJ_gamma,ADJ_lw,ADJ_zcld
  10    !          ADJ_theta (somtime theta is a variable)
  11    ! Input  : ADJ_tbup,ADJ_tbdn,ADJ_tauatm
  12    ! Output mean fields : tbup,tbdn,tauatm
  13    !-----------------------------------------------------------------
  14
  15    integer, intent(in)    :: ifreq
  16    real,    intent(in)    :: theta,p0,wv,hwv,ta,gamma,lw,zcld
  17    real,    intent(inout) :: ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta, ADJ_gamma,ADJ_lw,ADJ_zcld,ADJ_theta
  18    real,    intent(inout) :: ADJ_tbup,ADJ_tbdn,ADJ_tauatm
  19    real,    intent(out)   :: tbup,tbdn,tauatm
  20
  21    real :: tbdn_save
  22
  23    real :: mu,hdn,hup,hdninf,hupinf,ADJ_mu
  24
  25    real :: b1(4),b2(4),b3(4)
  26    real :: c(4),d1(4),d2(4),d3(4),zeta(4),kw0(4),kw1(4),kw2(4),kw3(4)
  27    real :: tau,tau1,tau2,taucld
  28    real :: tcld,tc,em,em1
  29    real :: sigv,sigo,sig,sig1,sigcld
  30    real :: teff1dn,teff1up,teffdn,teffup
  31    real :: tbcld,tbclrdn,tbclrup,tb1dn,tb1up,tb2dn,tb2up
  32    real :: otbar,tc2,tc3,hv,ho,alph
  33    real :: ADJ_sigv,ADJ_otbar,ADJ_sigo,ADJ_tcld,ADJ_tc,ADJ_tc2,ADJ_tc3
  34    real :: ADJ_sigcld,ADJ_taucld,ADJ_tbcld,ADJ_hv,ADJ_ho
  35    real :: ADJ_hdn,ADJ_hup,ADJ_hdninf,ADJ_sig,ADJ_sig1,ADJ_tau,ADJ_tau1
  36    real :: ADJ_tau2,ADJ_em1,ADJ_teff1dn,ADJ_hupinf,ADJ_em,ADJ_teff1up
  37    real :: ADJ_teffdn,ADJ_teffup,ADJ_tbclrdn,ADJ_tbclrup,ADJ_tb1dn,ADJ_tb1up
  38    real :: ADJ_tb2dn,ADJ_tb2up,ADJ_alph
  39
  40    data b1/-.46847e-1,-.57752e-1,-.18885,.10990/
  41    data b2/.26640e-4,.31662e-4,.9832e-4,.60531e-4/
  42    data b3/.87560e+1,.10961e+2,.36678e+2,-.37578e+2/
  43    data c/ .9207,   1.208,     .8253,     .8203/
  44    data zeta/4.2,4.2,4.2,2.9/
  45    data d1/-.35908e+1,-.38921e+1,-.43072e+1,-.17020e+0/
  46    data d2/ .29797e-1, .31054e-1, .32801e-1, .13610e-1/
  47    data d3/-.23174e-1,-.23543e-1,-.24101e-1,-.15776e+0/
  48    data kw0/ .786e-1, .103,    .267,    .988/
  49    data kw1/-.230e-2,-.296e-2,-.673e-2,-.107e-1/
  50    data kw2/ .448e-4, .557e-4, .975e-4,-.535e-4/
  51    data kw3/-.464e-6,-.558e-6,-.724e-6, .115e-5/
  52
  53    if (trace_use) call da_trace_entry("da_tbatmos_adj")
  54
  55    mu=0.0;hdn=0.0;hup=0.0;hdninf=0.0;hupinf=0.0;ADJ_mu=0.0
  56    tcld=0.0;tc=0.0;em=0.0;em1=0.0
  57    sigv=0.0;sigo=0.0;sig=0.0;sig1=0.0;sigcld=0.0
  58    teff1dn=0.0;teff1up=0.0;teffdn=0.0;teffup=0.0
  59    tbcld=0.0;tbclrdn=0.0;tbclrup=0.0;tb1dn=0.0;tb1up=0.0;tb2dn=0.0;tb2up=0.0
  60    otbar=0.0;tc2=0.0;tc3=0.0;hv=0.0;ho=0.0;alph=0.0
  61    ADJ_sigv=0.0;ADJ_otbar=0.0;ADJ_sigo=0.0;ADJ_tcld=0.0;
  62    ADJ_tc=0.0;ADJ_tc2=0.0;ADJ_tc3=0.0
  63    ADJ_sigcld=0.0;ADJ_taucld=0.0;ADJ_tbcld=0.0;ADJ_hv=0.0;ADJ_ho=0.0
  64    ADJ_hdn=0.0;ADJ_hup=0.0;ADJ_hdninf=0.0;ADJ_sig=0.0;ADJ_sig1=0.0
  65    ADJ_tau=0.0;ADJ_tau1=0.0
  66    ADJ_tau2=0.0;ADJ_em1=0.0;ADJ_teff1dn=0.0;ADJ_hupinf=0.0;ADJ_em=0.0
  67    ADJ_teff1up=0.0;ADJ_teffdn=0;ADJ_teffup=0.0;ADJ_tbclrdn=0.0
  68    ADJ_tbclrup=0.0;ADJ_tb1dn=0.0;ADJ_tb1up=0.0
  69    ADJ_tb2dn=0.0;ADJ_tb2up=0.0;ADJ_alph=0.0
  70    tau=0.0;tau1=0.0;tau2=0.0;taucld=0.0
  71    tcld=0.0;tc=0.0;em=0.0;em1=0.0
  72    sigv=0.0;sigo=0.0;sig=0.0;sig1=0.0;sigcld=0.0
  73    teff1dn=0.0;teff1up=0.0;teffdn=0.0;teffup=0.0
  74
  75
  76    ! mu = secant(theta)
  77    ! somtime theta is a variable
  78
  79    mu     = 1.0/cos(theta*0.0174533)
  80
  81    ! get water vapor optical depth
  82
  83    call cal_sigma_v(ifreq,p0,wv,hwv,ta,gamma,sigv)
  84
  85    ! otbar = one over "mean" temperature
  86
  87    otbar =   1.0/(ta - gamma*zeta(ifreq))
  88
  89    ! sigo = dry air optical depth
  90
  91    sigo = b1(ifreq) + b2(ifreq)*    p0  + b3(ifreq)*    otbar
  92
  93    ! cloud parameters
  94
  95    tcld   =     ta - gamma*zcld
  96    tc     =     tcld - t_kelvin
  97    tc2    =     tc*tc
  98    tc3    =     tc2*tc
  99    sigcld =  ( kw0(ifreq) + tc*kw1(ifreq) + tc2*kw2(ifreq) +  &
 100                tc3*kw3(ifreq) )*lw
 101    taucld =   exp(-mu*sigcld)
 102    tbcld  =   (1.0 - taucld)*tcld
 103
 104    ! hv, ho = effective absorber scale heights for vapor, dry air
 105
 106    hv = c(ifreq)*   hwv
 107    ho = d1(ifreq) + d2(ifreq)* ta + d3(ifreq)* gamma
 108
 109    ! get effective emission heights for layer 1 and total atmosphere
 110
 111
 112    call effht(ho,hv,sigo,sigv,mu,zcld,hdn,hup, hdninf,hupinf)
 113
 114    ! atmospheric transmittances in layer one and two, and combined
 115
 116    sig =     sigo +     sigv
 117    sig1 = sigo*(1.0-exp(-zcld/ho)) + sigv*(1.0-exp(-zcld/hv))
 118    tau  =  exp(-mu*sig)
 119    tau1 =  exp(-mu*sig1)
 120    tau2 =  tau/tau1
 121
 122    ! atmospheric "emissivity"
 123
 124    em1  =   1.0 - tau1
 125    em   =   1.0 - tau
 126
 127    ! downwelling and upwelling brightness temperature for each layer
 128
 129    teff1dn =     ta - gamma*hdn
 130    teff1up =     ta - gamma*hup
 131    teffdn =     ta - gamma*hdninf
 132    teffup =     ta - gamma*hupinf
 133    tbclrdn = teffdn*em
 134    tbclrup = teffup*em
 135
 136    tb1dn = em1*teff1dn
 137    tb1up = em1*teff1up
 138    tb2dn = (tbclrdn - tb1dn)/tau1
 139    tb2up =      tbclrup - tau2*tb1up
 140
 141    ! total downwelling and upwelling brightness temperature and transmittance
 142
 143    tbdn  =     tb1dn + tau1*(tbcld + taucld*tb2dn)
 144    tbup  =     tb2up + tau2*(tbcld + taucld*tb1up)
 145    tauatm = tau*taucld
 146
 147    ! the following lines apply an ad hoc correction to improve fit
 148    ! at large angles and/or high gaseous opacities
 149    !  (downwelling brightness temperatures only)
 150
 151    alph = (0.636619*atan(mu*sig))**2
 152    tbdn_save = tbdn
 153    tbdn = (1.0-alph)*tbdn + em*alph*ta
 154
 155
 156    ! start
 157
 158    tbdn = tbdn_save
 159
 160    ADJ_alph = - ADJ_tbdn*tbdn
 161    ADJ_em   = ADJ_tbdn*alph*ta
 162    ADJ_alph = em*ADJ_tbdn*ta + ADJ_alph
 163    ADJ_ta   = em*alph*ADJ_tbdn + ADJ_ta
 164    ADJ_tbdn = (1.0-alph)*ADJ_tbdn
 165
 166    if (abs(sig) .gt. 0.0) then
 167       ADJ_mu  = 2.0*0.636619*0.636619*ADJ_alph*sig*atan(mu*sig)/(1.0+mu*mu*sig*sig)
 168       ADJ_sig = 2.0*0.636619*0.636619*mu*ADJ_alph*atan(mu*sig)/(1.0+mu*mu*sig*sig)
 169    else
 170       ADJ_mu  = 0.0
 171       ADJ_sig = 0.0
 172    end if
 173
 174    ADJ_tau    = ADJ_tauatm*taucld
 175    ADJ_taucld = tau*ADJ_tauatm
 176    ADJ_tb2up  = ADJ_tbup
 177    ADJ_tau2   = ADJ_tbup*(tbcld + taucld*tb1up)
 178    ADJ_tbcld  = tau2*ADJ_tbup
 179    ADJ_taucld = tau2*ADJ_tbup*tb1up  + ADJ_taucld
 180    ADJ_tb1up  = tau2*taucld*ADJ_tbup
 181    ADJ_tb1dn = ADJ_tbdn
 182    ADJ_tau1  = ADJ_tbdn*(tbcld + taucld*tb2dn)
 183    ADJ_tbcld = tau1*ADJ_tbdn   + ADJ_tbcld
 184    ADJ_taucld = tau1*ADJ_tbdn*tb2dn  + ADJ_taucld
 185    ADJ_tb2dn  = tau1*taucld*ADJ_tbdn
 186
 187    ADJ_tbclrup =   ADJ_tb2up
 188    ADJ_tau2    = - ADJ_tb2up*tb1up + ADJ_tau2
 189    ADJ_tb1up   = - tau2*ADJ_tb2up  + ADJ_tb1up
 190
 191    ADJ_tbclrdn =   ADJ_tb2dn/tau1
 192    ADJ_tb1dn   = - ADJ_tb2dn/tau1  + ADJ_tb1dn
 193    ADJ_tau1    = - tb2dn*ADJ_tb2dn/tau1 + ADJ_tau1
 194
 195    ADJ_em1     = ADJ_tb1up*teff1up
 196    ADJ_teff1up = em1*ADJ_tb1up
 197
 198    ADJ_em1     = ADJ_tb1dn*teff1dn + ADJ_em1
 199    ADJ_teff1dn = em1*ADJ_tb1dn
 200
 201    ADJ_teffup  = ADJ_tbclrup*em
 202    ADJ_em      = teffup*ADJ_tbclrup + ADJ_em
 203
 204    ADJ_teffdn  = ADJ_tbclrdn*em
 205    ADJ_em      = teffdn*ADJ_tbclrdn + ADJ_em
 206
 207    ADJ_ta      =   ADJ_teffup + ADJ_ta
 208    ADJ_gamma   = - ADJ_teffup*hupinf + ADJ_gamma
 209    ADJ_hupinf  = - gamma*ADJ_teffup
 210
 211    ADJ_ta      =   ADJ_teffdn + ADJ_ta
 212    ADJ_gamma   = - ADJ_teffdn*hdninf + ADJ_gamma
 213    ADJ_hdninf  = - gamma*ADJ_teffdn
 214
 215    ADJ_ta    =   ADJ_teff1up + ADJ_ta
 216    ADJ_gamma = - ADJ_teff1up*hup + ADJ_gamma
 217    ADJ_hup   = - gamma*ADJ_teff1up
 218
 219    ADJ_ta    =   ADJ_teff1dn + ADJ_ta
 220    ADJ_gamma = - ADJ_teff1dn*hdn + ADJ_gamma
 221    ADJ_hdn   = - gamma*ADJ_teff1dn
 222
 223    ADJ_tau  = - ADJ_em + ADJ_tau
 224
 225    ADJ_tau1 = - ADJ_em1 + ADJ_tau1
 226
 227    ADJ_tau  =   ADJ_tau2/tau1 + ADJ_tau
 228    ADJ_tau1 = - tau2*ADJ_tau2/tau1 + ADJ_tau1
 229
 230    ADJ_sig1 = - mu*ADJ_tau1*tau1
 231    ADJ_mu   = - ADJ_tau1*sig1*tau1 + ADJ_mu
 232
 233    ADJ_mu   = - ADJ_tau*sig*tau + ADJ_mu
 234    ADJ_sig  = - mu*ADJ_tau*tau + ADJ_sig
 235
 236    ADJ_sigo = ADJ_sig1*(1.0-exp(-zcld/ho))
 237    ADJ_sigv = ADJ_sig1*(1.0-exp(-zcld/hv))
 238    ADJ_zcld = sigo*ADJ_sig1/ho*exp(-zcld/ho) + ADJ_zcld
 239    ADJ_ho   = - sigo*zcld*ADJ_sig1/(ho*ho)*exp(-zcld/ho)
 240    ADJ_zcld = sigv*ADJ_sig1/hv*exp(-zcld/hv) + ADJ_zcld
 241    ADJ_hv   = - sigv*zcld*ADJ_sig1/(hv*hv)*exp(-zcld/hv)
 242
 243    ADJ_sigo = ADJ_sig + ADJ_sigo
 244    ADJ_sigv = ADJ_sig + ADJ_sigv
 245
 246    call da_effht_adj(ho,hv,sigo,sigv,mu,zcld,hdn,hup,        &
 247                   hdninf,hupinf,                          &
 248                   ADJ_ho,ADJ_hv,ADJ_sigo,ADJ_sigv,ADJ_mu, &
 249                   ADJ_zcld,ADJ_hdn,ADJ_hup,ADJ_hdninf,    &
 250                   ADJ_hupinf                              )
 251
 252    ADJ_ta    = d2(ifreq)*ADJ_ho + ADJ_ta
 253    ADJ_gamma = d3(ifreq)*ADJ_ho + ADJ_gamma
 254    ADJ_hwv   = c(ifreq)*ADJ_hv  + ADJ_hwv
 255
 256    ADJ_taucld = - ADJ_tbcld*tcld + ADJ_taucld
 257    ADJ_tcld = (1.0 - taucld)*ADJ_tbcld
 258
 259    ADJ_mu     = - ADJ_taucld*sigcld*taucld + ADJ_mu
 260    ADJ_sigcld = - mu*ADJ_taucld*taucld
 261
 262    ADJ_tc  = ADJ_sigcld*kw1(ifreq)*lw
 263    ADJ_tc2 = ADJ_sigcld*kw2(ifreq)*lw
 264    ADJ_tc3 = ADJ_sigcld*kw3(ifreq)*lw
 265    ADJ_lw  = (kw0(ifreq)+tc*kw1(ifreq)+tc2*kw2(ifreq)+tc3*kw3(ifreq)) &
 266              *ADJ_sigcld + ADJ_lw
 267
 268    ADJ_tc2  = ADJ_tc3*tc + ADJ_tc2
 269    ADJ_tc   = tc2*ADJ_tc3 + ADJ_tc
 270    ADJ_tc   = 2.0*tc*ADJ_tc2 + ADJ_tc
 271
 272    ADJ_tcld = ADJ_tc + ADJ_tcld
 273
 274    ADJ_ta    =   ADJ_tcld + ADJ_ta
 275    ADJ_gamma = - ADJ_tcld*zcld + ADJ_gamma
 276    ADJ_zcld  = - gamma*ADJ_tcld + ADJ_zcld
 277
 278    ADJ_p0    = b2(ifreq)*ADJ_sigo + ADJ_p0
 279    ADJ_otbar = b3(ifreq)*ADJ_sigo
 280
 281    ADJ_ta    = - otbar*otbar*ADJ_otbar + ADJ_ta
 282    ADJ_gamma =   otbar*otbar*ADJ_otbar*zeta(ifreq) + ADJ_gamma
 283
 284    call da_sigma_v_adj(ifreq,p0,wv,hwv,ta,gamma,sigv,   &
 285                         ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta,    &
 286                         ADJ_gamma,ADJ_sigv               )
 287
 288    ADJ_theta = mu*mu*0.0174533*ADJ_mu*sin(theta*0.0174533) + ADJ_theta
 289
 290    if (trace_use) call da_trace_exit("da_tbatmos_adj")
 291
 292 end subroutine da_tbatmos_adj
 293
 294