subroutine get_profile_data_cfsr(debug, rlat, met_year, met_month, * met_day, land, theta, sat, psfc, pmsl, tp, w, o3mr, * tozone, shifted_FM_opt, z, isp, taup, sfc_emis, ptrp, * ltrp, ttpp, rclr_s, rclr_s2) c!F77------------------------------------------------------------------- c!Description: Get profile and related surface data for current pixel c c!Input parameters: c c debug Debug write flag c rlat Latitude c met_year Year c met_month Month c met_day Day of year c land Land/sea flag (0=water, 1=land) c theta Viewing zenith angle c sat NOAA spacecraft c psfc Surface pressure c pmsl Pressure at mean sea level c tp 101-level temperature profile c w 101-level water vapor mixing ratio profile c shifted_FM_opt shifted/unshifted forward model option c 1=unshifted (original), 2=shifted FM and Planck coefficients c o3mr 6-level O3 profile from ancillary data c tozone Total column O3 from ancillary data c c!Output parameters: c z 101-level geopotential height profile c isp Level (1-101) of surface pressure c taup 101-level Transmittance profiles c sfc_emis Estimate of surface emissivity c ptrp Atm. pressure at tropopause c ltrp Level (1-101) of tropopause c ttpp Brightness temperature profile at 11 microns c assuming "black" surface c rclr_s 11 micron radiance profile corresponding to 'ttpp' c rclr_s2 12 micron radiance profile corresponding to 'ttpp' c c!Revision history: c 07/10 R. Frey Original version c 06/11 R. Frey Added 'shifted_FM_opt' c c c!Calls: c clozo101 c height c adjo3 c rf_hirstran_101 c assign_eco_emis c getiremis c function rf_hirsbright c function fmrad_emis c c!End------------------------------------------------------------------- implicit none save integer plev parameter (plev = 101) integer ntbct parameter (ntbct = 8) integer kwc parameter (kwc = 5) integer nb_wavelen parameter (nb_wavelen = 7) integer nb, nt parameter (nb=20,nt=2) real cwn(nb), fk1(nb), fk2(nb), tc(nt,nb) real cwn_shifted(nb), fk1_shifted(nb), fk2_shifted(nb), tc_shifted(nt,nb) common/hirpfc/cwn, fk1, fk2, tc common/hirpfc_shifted/cwn_shifted, fk1_shifted, fk2_shifted, tc_shifted c Scalar arguments. character*6 sat integer debug, met_year, met_month, met_day, isp, ltrp, land, * shifted_FM_opt real rlat, theta, psfc, ptrp, pmsl, tozone c Array arguments. real z(plev), tp(plev), w(plev), taup(plev,ntbct), * sfc_emis(ntbct), ttpp(plev), rclr_s(plev), rclr_s2(plev), * o3mr(6) c Local scalars. integer ii, jj, ll, k, iisp, nl, lmin, kban, jdet, lsf, * imslp, init real zs, emis_out, rho_out, tss, tmin, emiswc, emis12, * ppsfc c Local arrays. real pp(plev), ozone(plev), adj_ozone(plev), * emis(nb_wavelen), freqemis(nb_wavelen), rho(nb_wavelen), * freq(ntbct) integer mbnds(ntbct), std_doy(11), leap_doy(11) c External functions. real fmrad_emis, rf_hirsbright external fmrad_emis, rf_hirsbright 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 ... Initialize ozone profile. data ozone/plev*0.0/ c Dectector number for transmittance calculations (0=averaged). data jdet /0/ c Day of year at end of first 11 months: standard and leap year. data std_doy /31,59,90,120,151,181,212,243,273,304,334/ data leap_doy /31,60,91,121,152,182,213,244,274,305,335/ c----------------------------------------------------------------------- c Perform one-time initializations. if(init .eq. 0) then c Choose unshifted (original=1) or shifted (=2) FM and Planck coefficients. c Store central wave numbers of desired channels (see 'mbnds'). if(shifted_FM_opt .eq. 1) then do k = 1, ntbct freq(k) = cwn(mbnds(k)) enddo else if(shifted_FM_opt .eq. 2) then do k = 1, ntbct freq(k) = cwn_shifted(mbnds(k)) enddo end if init = 1 end if c----------------------------------------------------------------------- c Initialize geopotential height profile. do k = 1, plev z(k) = 0.0 enddo c---------------------------------------------------------------------- c Get climatological ozone profile. call clozo101(rlat, met_month, ozone) c---------------------------------------------------------------------- c Adjust to ozone profile according to GDAS data. call adjo3(debug, ozone, o3mr, tozone, adj_ozone) c---------------------------------------------------------------------- c Find level of surface pressure. do ll = 1,plev if(psfc .le. pp(ll)) then isp = ll go to 100 end if enddo 100 continue c---------------------------------------------------------------------- c Find level of mean sea level pressure. do ll = 1,plev if(pmsl .le. pp(ll)) then imslp = ll go to 200 end if enddo 200 continue c---------------------------------------------------------------------- c Get geopotential height profile (km). zs = 0.0 nl = imslp call height(pp,tp,w,zs,nl,z) c---------------------------------------------------------------------- c Calculate transmittance profiles (101 level fast model) do k = 1, ntbct kban = mbnds(k) call rf_hirstran_101(debug, met_year, met_day, tp, w, adj_ozone, & theta, sat, kban, shifted_FM_opt, & taup(1,k)) enddo c---------------------------------------------------------------------- c Write debug information. if(debug .gt. 1) then write(*,'(''Atmospheric transmittance profiles:'')') write(*,'(''Profile Date: '',3i10,a10)') met_year, met_month, & met_day, sat write(*,'(''Atm level, pressure, temperature, '', & ''water vapor, ozone, adj_ozone, transmittances (4,5,6,7,8,10,11,12)'')') do k=1,plev write(*,'(i4,14f10.4)') k,pp(k),tp(k),w(k),ozone(k),adj_ozone(k), & (taup(k,jj),jj=1,ntbct) enddo write(*,'(1x)') end if c---------------------------------------------------------------------- c Get surface emissivity. call assign_eco_emis(land,emis,rho,freqemis) do k = 1, ntbct call getiremis(nb_wavelen,freqemis,freq(k),emis,rho,emis_out, & rho_out) sfc_emis(k) = emis_out enddo if(debug .gt. 1) then write(*,'(''Surface emissivity calculations'')') write(*,'(''Land/sea tag: '',i10)') land write(*,'(''Reference frequencies: '',7f10.3)') freqemis write(*,'(''Reference emissivities: '',7f10.3)') emis write(*,'(''Reference rho: '',7f10.3)') rho write(*,'(''Instrument freqencies: '',10f10.3)') freq write(*,'(''Output emissivities: '',10f10.3)') sfc_emis write(*,'(1x)') end if c---------------------------------------------------------------------- c Find level of coldest temperature between 100 mb and the surface. tmin = 350.0 do ll = 44,isp if(tp(ll) .le. tmin) then tmin = tp(ll) lmin = ll end if enddo c---------------------------------------------------------------------- c Check tropopause height. Look for inflection point in temperature profile c if found do not change trop. level. c In isothermal conditions, choose the bottom level for trop. intead of top level. c If temperature increases monotonically below 100 mb, do not change trop level. ltrp = -99 do ll = lmin, 71 if(tp(ll-1) .ge. tp(ll) .and. tp(ll+1) .gt. tp(ll)) then ltrp = ll go to 300 end if enddo 300 continue if(ltrp .eq. -99) ltrp = lmin ptrp = pp(ltrp) c---------------------------------------------------------------------- c Get IRW radiance profile. do k = ltrp,isp if(k .eq. isp) then emiswc = sfc_emis(kwc) emis12 = sfc_emis(6) ppsfc = psfc else emiswc = 1.0 emis12 = 1.0 ppsfc = pp(k) end if tss = tp(k) iisp = k rclr_s(k) = fmrad_emis(tss,emiswc,ppsfc,pp,tp,taup(1,kwc), & mbnds(kwc),iisp,shifted_FM_opt) TTPP(k) = rf_hirsbright(rclr_s(k),mbnds(kwc), shifted_FM_opt) rclr_s2(k) = fmrad_emis(tss,emis12,ppsfc,pp,tp,taup(1,6), & mbnds(6),iisp,shifted_FM_opt) c write(*,'(''ttpp: '', i5,4f12.6)') k, pp(k), rclr_s(k), c & rclr_s2(k), ttpp(k) enddo if(debug .gt. 1) then write(*,'(1x)') write(*,'(''Calculated brightness temperature and radiance '', & ''profiles from trop to sfc:'')') write(*,'(''Atm. level, Atm. pressure, Tbb, 11 um radiance, '', & ''12 um radiance'')') do k = ltrp, isp write(*,'(i5, 4f12.4)') k, pp(k), ttpp(k), rclr_s(k), rclr_s2(k) enddo write(*,'(1x)') end if c----------------------------------------------------------------------- return end