subroutine CO2_retrieval(debug, landsea, beta, wtrcld_frc, delr, ra, ltrp, * iw1, tp, taup, rwarm, ecawin, shifted_FM_opt, * pct, eca, tct, ip, lco2, co2_flag) c!F77------------------------------------------------------------------- c!Description: Calculate CO2-slicing cloud top pressure, effective cloud c amount, and cloud top temperature. c c!Input parameters: c debug Debug write flag c landsea Land/sea flag c beta Beta ratio (11/12 microns) c delr LHS of CO2-slicing equation c wtrcld_frc Collocated PATMOS-X (AVHRR) water cloud fraction c in the current HIRS FOV c ra RHS of CO2-slicing equation c ltrp Level of tropopause (1-101) c iw1 Lowest level for CO2-slicing result c tp Temperature profile c taup Transmittance profiles c rwarm Calculated clear-sky radiances c ecawin Cloud amount from IRW retrieval c shifted_FM_opt shifted/unshifted forward model option c 1=unshifted (original), 2=shifted FM and Planck c coefficients c c!Output parameters: c pct Cloud top pressure c eca Effective cloud amount c tct Cloud top temperature c ip Index of band pair chosen for final result c lco2 Level of cloud top pressure retrieval (1-101) c co2_flag Indicates successful CO2-slicing retrieval (=2) c c!Revision history: c 07/10 R. Frey Original version c 06/11 R. Frey Added 'shifted_FM_opt' c 12/13 R. Frey Added 'wtrcld_frc' c c!Calls: c function fmrad_emis c c!End------------------------------------------------------------------- implicit none save integer plev parameter (plev = 101) integer nbct parameter (nbct = 5) integer ntbct parameter (ntbct = 8) integer nsct parameter (nsct = 4) integer kwc parameter (kwc = 5) c Scalar arguments. integer debug, landsea, ltrp, iw1, ip, lco2, co2_flag, * shifted_FM_opt real beta, ecawin, pct, eca, tct, wtrcld_frc c Array arguments. real delr(nbct), ra(plev,nbct), tp(plev), taup(plev,ntbct), * rwarm(nbct) c Local scalars. integer id, k1, k2, ll, ipco2, ngch, k, ii real beta_threshold, fmsav, fm1, fm2, * fm, bot, top, pfco2, emisrw, bad_value logical ok, neg, start c Local arrays. integer krto(nsct), kch(nsct,2), lev(nsct), mbnds(ntbct) real rmin(nbct), rwcld(nsct), amo(nsct), ratio(nsct), pp(plev), * em_adj(nsct) c External functions. real fmrad_emis external fmrad_emis c----------------------------------------------------------------------- c Define 101 pressure levels. data pp / 0.0050, 0.0161, 0.0384, 0.0769, 0.1370, + 0.2244, 0.3454, 0.5064, 0.7140, 0.9753, 1.2972, + 1.6872, 2.1526, 2.7009, 3.3398, 4.0770, 4.9204, + 5.8776, 6.9567, 8.1655, 9.5119, 11.0038, 12.6492, + 14.4559, 16.4318, 18.5847, 20.9224, 23.4526, 26.1829, + 29.1210, 32.2744, 35.6505, 39.2566, 43.1001, 47.1882, + 51.5278, 56.1260, 60.9895, 66.1253, 71.5398, 77.2396, + 83.2310, 89.5204, 96.1138, 103.0172, 110.2366, 117.7775, + 125.6456, 133.8462, 142.3848, 151.2664, 160.4959, 170.0784, + 180.0183, 190.3203, 200.9887, 212.0277, 223.4415, 235.2338, + 247.4085, 259.9691, 272.9191, 286.2617, 300.0000, 314.1369, + 328.6753, 343.6176, 358.9665, 374.7241, 390.8926, 407.4738, + 424.4698, 441.8819, 459.7118, 477.9607, 496.6298, 515.7200, + 535.2322, 555.1669, 575.5248, 596.3062, 617.5112, 639.1398, + 661.1920, 683.6673, 706.5654, 729.8857, 753.6275, 777.7897, + 802.3714, 827.3713, 852.7880, 878.6201, 904.8659, 931.5236, + 958.5911, 986.0666, 1013.9476, 1042.2319, 1070.9170, 1100.0000/ c Order of HIRS IR band numbers used data mbnds/4,5,6,7,8,10,11,12/ c Beta ratio threshold. c data beta_threshold /0.95/ data beta_threshold /0.80/ c Define noise level (in radiance units -- mw/m2/ster/cm-1) for each c CO2 channel. c data rmin /-1.0, -1.0, -1.0, -1.0, -1.0/ c C6 Aqua values for testing. c data rmin /-1.25, -1.0, -4.0, -4.0, -0.5/ data rmin /-0.5, -0.5, -0.5, -0.5, -0.5/ c Define indices of channel combinations used in the co2-slicing c method. data kch /1,2,2,3,2,3,4,4/ c Ice cloud emissivity adjustments for each channel combination. data em_adj / 1.0, 1.0, 1.0, 1.0 / c Black cloud emissivity data emisrw /1.0/ c Fill value. data bad_value /-999.99/ c----------------------------------------------------------------------- c Initializations. pct = -99.9 eca = -99.9 tct = -99.9 ip = -99 lco2 = -99 ngch = 0 co2_flag = 1 do id = 1, nsct krto(id) = 0 lev(id) = 0 enddo c----------------------------------------------------------------------- c Compute cloud top pressure for each band pair when possible. c Check for water clouds over water surfaces. c if(landsea .eq. 0) then c if(beta .lt. beta_threshold) go to 3000 c endif do id = 1, nsct krto(id) = 0 fmsav = 1000.0 k1 = kch(id,1) k2 = kch(id,2) c ... Check if values of cold minus warm are within instrument noise. if(delr(k1) .le. rmin(k1)) then if(delr(k2) .le. rmin(k2)) then c ... Find intersection between LHS and RHS of co2-slicing c ... equation. Apply emissivity correction. ok = .false. do ll = ltrp, iw1 fm1 = delr(k1) / delr(k2) fm2 = (ra(ll,k1) * em_adj(id)) / ra(ll,k2) fm = fm1 - fm2 if(fm .lt. 0.0) then neg = .true. else neg = .false. end if if(ll .eq. ltrp) start = neg if(neg .neqv. start) then ok = .true. lev(id) = ll - 1 fmsav = abs(fm) go to 2000 end if if(debug .gt. 2) then write(*,'(4i5,3f12.5,2l5)') id,k1,k2,ll,fm1,fm2,fm, & neg,start end if enddo 2000 continue if(fmsav .lt. 1000.0) then krto(id) = 1 end if if( (.not. ok) ) krto(id) = 0 if( (.not. ok) .and. neg) then krto(id) = 1 lev(id) = iw1 end if end if end if enddo c----------------------------------------------------------------------- if(debug .gt. 1) then write(*,'(1x)') write(*,'(''Initial CO2-slicing cloud levels'')') write(*,'(''HIRS channel combinations:'')') write(*,'(9x,8(i8,''/'',i1))') (mbnds(kch(ii,1)), & mbnds(kch(ii,2)),ii=1,nsct) write(*,'(''atm lev '',8i10)') (lev(ii),ii=1,nsct) write(*,'(1x)') end if c----------------------------------------------------------------------- c ... Compute effective cloud amounts for every channel combination c ... for which there was a successful cloud height retrieval. do id = 1, nsct rwcld(id) = bad_value amo(id) = bad_value if(krto(id) .ne. 0) then ll = lev(id) rwcld(id) = fmrad_emis(tp(ll),emisrw,pp(ll),pp,tp, & taup(1,kwc),mbnds(kwc),ll,shifted_FM_opt) bot = rwcld(id) - rwarm(kwc) top = delr(kwc) if(abs(bot) .gt. 0.1) then ratio(id) = top / bot if(ratio(id) .le. 1.0 .and. ratio(id) .ge. 0.01) then amo(id) = ratio(id) * ecawin else krto(id) = 0 end if else krto(id) = 0 end if end if enddo c---------------------------------------------------------------------- c Write debug information about effective emissivity calculations if(debug .gt. 1) then write(*,'(''Initial CO2-slicing effective emissvities'')') write(*,'(''HIRS channel combinations:'')') write(*,'(9x,8(i8,''/'',i1))') (mbnds(kch(ii,1)), & mbnds(kch(ii,2)),ii=1,nsct) write(*,'(''cld rad '',8f10.3)') & (rwcld(ii),ii=1,nsct) write(*,'(''ratios '',8f10.3)') & (ratio(ii),ii=1,nsct) write(*,'(''eff emiss '',8f10.3)') & (amo(ii),ii=1,nsct) write(*,'(1x)') end if c---------------------------------------------------------------------- c ... Use "top-down" method for choosing a solution. if(krto(4) .ne. 0) then if(pp(lev(4)) .lt. 650.0) then ipco2 = 4 lco2 = lev(4) ngch = 1 end if end if if(krto(2) .ne. 0) then if(pp(lev(2)) .lt. 550.0) then ipco2 = 2 lco2 = lev(2) ngch = 1 end if end if if(krto(1) .ne. 0) then if(pp(lev(1)) .lt. 450.0) then ipco2 = 1 lco2 = lev(1) ngch = 1 end if end if c---------------------------------------------------------------------- c Write debug information about final CTP selection. if(debug .gt. 1) then write(*,'(''CO2-slicing cloud top pressure selection'')') write(*,'(''HIRS channel combinations:'')') write(*,'(9x,8(i8,''/'',i1))') (mbnds(kch(ii,1)), & mbnds(kch(ii,2)),ii=1,nsct) write(*,'(''heights '',8f10.3)') (pp(lev(ii)),ii=1,nsct) write(*,'(''eff emiss '',8f10.3)') (amo(ii),ii=1,nsct) write(*,'(1x)') end if c---------------------------------------------------------------------- c Check if there were any valid co2-slicing cloud heights. c Perform additional check on AVHRR water cloud fraction. if(ngch .ne. 0) then c Best CO2-slicing cloud top pressure estimate. pfco2 = pp(lco2) c Check for possible mid-level water phase cloud. if(landsea .ne. 0 .or. (wtrcld_frc .lt. 0.75 .or. pfco2 .lt. 440.0)) then c Prepare final cloud height for output and set appropritate c flags. c Round cloud pressure to nearest 5 mb. pct = (nint(pfco2 / 5.0)) * 5.0 eca = amo(ipco2) ip = ipco2 tct = tp(lco2) co2_flag = 2 end if end if c----------------------------------------------------------------------- 3000 continue return end