subroutine reflec(irefl,iozon,iofset,ilat,xlat,isnow,grref,clref, 1 sza12,pcloud,pteran,estozn,xnvalm,xnvalp,ramcor,fteran,frstrf, 2 lprint,clfrac,ref,pref,iloprf,ihiprf,pcfrac,dnd331,dndrm, 3 dndrp) c c*********************************************************************** c c reflec c c august, 2001 by charlie wellemeyer of raytheon corporation c c purpose c compute cloud fraction and reflectivity for 331 nm wavelength c c method c establish assumed ground and cloud reflectivities. c select interpolation indices for table lookup. c perform table interpolations in pressure and ozone, c then calculate cloud fraction. check for the prsence c of snow or sunglint, for which clear sky is assumed. c calculate reflectivity by inverting the surface c radiance formula for clear or completly cloudy, and c use the partial cloud model for partly clouded scenes. c c variables c name type i/o description c ---- ---- --- ----------- c arguments c irefl i*4 i index for reflectivity wavelength c iozon i*4 i index for ozone wavelength c iofset i*4 i offset pointer. c ilat i*4 i latitude index (1 - low 2 - mid 3 - high) c xlat i*4 i latitude c isnow i*4 i snow/ice indicator c grref r*4 i/o ground reflectivity c clref r*4 i/o cloud reflectivity c sza12 r*4 i solar zenith angles for first 12 wavelengths c pcloud r*4 i cloud top pressure c pteran r*4 i terrain pressure c estozn r*4 i current total ozone estimate c xnvalm r*4 i measured n-values from the monochromator c xnvalp r*4 i measured n-values from the photometer c ramcor r*4 i raman scattering corrections c fteran r*4 i interpolation fraction for terrain pressure c frstrf l*4 i first call indicator c lprint l*4 i if = .true., print c clfrac r*4 o cloud fraction c ref r*4 o reflectivity for wavelength irefl c pref r*4 o perturbed ref (needed for dN/dR calculation) c pcfrac r*4 o perturbed cloud frac(needed for dN/dR calculation) c refpo r*4 o phot refl at ozone position of mono c dnd331 r*4 o ozone sensitivity of 331 nm channel c dndrm r*4 o reflectivity sensitivity of 331 nm channel c c*********************************************************************** c c -- input parameters implicit none real tgrr,tcll,ezigr,ezicl real xlat, pcloud, pteran, estozn, ref, pref real grref, clref, xnvalm(12), xnvalp(12), sza12(12) real ramcor(9,6), nval integer iofset(12), isnow, irefl, iozon, ilat, iloprf, ihiprf logical frstrf, glint, lprint(30) c -- internal parameters real ezero, t, qgclc(4), qcclc(4), qsavlo, qsavhi, alb real den, fteran, oznind, glnfrc, omeglo, omeghi, ozfrac, 1 rmcrgr, rmcrcl, ezlggr, tloggr, sbgrnd, grad, 2 ezlgcl, tlogcl, sbcld, crad, r331, pden, pgref, pcref, 3 clfrco, refpo, dndrp, clfrcr, pxnvalm, palb real ezgr(2),ssgr(2),tgr(2),sbgr(2),ezcl(2),sscl(2), 1 tcl(2),sbcl(2),radgr(2),radcl(2),qgcalc(2),qccalc(2),sb integer iprfl, ioz c -- output parameters real clfrac, pgrref, pclref, pcfrac, dnd331, dndrm c -- common area and data statements include 'stndprof.com' c real fgtmp, stprf, profl, terroz c common /stndprof/fgtmp(11),stprf(11,21), c 1 profl(21),terroz(21) c c -- convert reflectivity channel n value to albedo c (adjusted for scene change to the ozone channel) c c alb = 10.**(-xnvalm(irefl)) nval = (xnvalm(irefl)+(xnvalp(iozon)-xnvalp(irefl))) alb = 10.**(-nval) c c -- initialize ground and cloud reflectivities c call prflec(pteran,pcloud,lprint,grref,clref) c if (lprint(11)) then write (6,1000) xlat, isnow, iofset(1), grref, clref, 1 irefl, pcloud, pteran, alb, ilat, sza12(irefl), 2 estozn, 3 fteran, frstrf, ramcor endif c c -- perturb measured reflectivity channel n value by 0.1% c -- all varaibles beginning with p are from perturbed n value c -- (needed for dN/dR calculation) c pxnvalm = 1.002*nval palb = 10.**(-pxnvalm) pgrref = grref pclref = clref c c -- determine profile indices for first call, otherwise c -- use profile indeces determined by subroutine oznot c if (frstrf) then call prfind(xlat,estozn,lprint,ilat,iloprf,ihiprf) endif iprfl = iloprf - 2 c c -- do two calculations for bracketting ozone profiles c do ioz=1,2 iprfl = iprfl + 1 c c -- compute i0, tr, and sb values c call iztrsb(irefl,xlat,iprfl,pteran,pcloud, 1 iofset(irefl),lprint,ezgr(ioz),ssgr(ioz),tgr(ioz), 2 sbgr(ioz),ezcl(ioz),sscl(ioz),tcl(ioz),sbcl(ioz)) c c -- store radiance for ozone sensitivity calculation c radgr(ioz) = ezgr(ioz)+grref*tgr(ioz)/(1.0-grref*sbgr(ioz)) radcl(ioz) = ezcl(ioz)+clref*tcl(ioz)/(1.0-clref*sbcl(ioz)) c enddo c c -- compute delta ozone including terrain height corrections c omeglo = profl(iloprf) - terroz(iloprf)*fteran omeghi = profl(ihiprf) - terroz(ihiprf)*fteran ozfrac = (estozn-omeglo)/(omeghi-omeglo) c c -- calculate raman corrections for ground and cloud c call getrng(irefl,pteran,pcloud,grref,clref,ramcor,lprint, 1 rmcrgr,rmcrcl) c c -- Interpolate log of table parameters to current ozone estimate c -- and convert to radiance for cloud fraction calculation c ezlggr=alog10(ezgr(1))+alog10(ezgr(2)/ezgr(1))*ozfrac tloggr=alog10(tgr(1))+alog10(tgr(2)/tgr(1))*ozfrac sbgrnd=sbgr(1)+(sbgr(2)-sbgr(1))*ozfrac grad=10.0**ezlggr+grref*10.0**tloggr/(1.0-grref*sbgrnd) grad = grad + grad * rmcrgr c ezlgcl=alog10(ezcl(1))+alog10(ezcl(2)/ezcl(1))*ozfrac tlogcl=alog10(tcl(1))+alog10(tcl(2)/tcl(1))*ozfrac sbcld=sbcl(1)+(sbcl(2)-sbcl(1))*ozfrac crad=10.0**ezlgcl+clref*10.0**tlogcl/(1.0-clref*sbcld) crad = crad + crad * rmcrcl c c -- Calculate cloud fraction c if (crad-grad.ne.0.) then clfrac = (alb-grad) / (crad-grad) pcfrac = (palb-grad) / (crad-grad) else clfrac = 0. pcfrac = 0. endif c c assume clear sky if snow is present (v7 used ge 5 total.f) c c if(isnow.eq.10) clfrac = 0.0 if(isnow.ge.5) clfrac = 0.0 c c -- calculate reflectivity using version 6 method for cloud free c if (clfrac.le.0.0) then ezigr=10.0**(ezlggr) ezigr=ezigr+ezlggr*rmcrgr den = alb - ezigr pden = palb - ezigr if (den.ne.0.) then tgrr=10.0**(tloggr) tgrr=tgrr+tgrr*rmcrgr ref = 1. / (tgrr/den + sbgrnd) else ref = 0. endif if (pden.ne.0.) then tgrr=10.0**(tloggr) tgrr=tgrr+tgrr*rmcrgr pref = 1. / (tgrr/pden + sbgrnd) else pref = 0. endif clfrac = 0.0 pcfrac = 0.0 grref = ref pgref = pref c if(frstrf) dnd331 = (-100.0*alog10(radgr(2)/radgr(1)))/ 1 (omeghi-omeglo) c c -- calculate reflectivity using version 6 method for cloudy case c else if (clfrac.ge.1.) then ezicl=10.0**(ezlgcl) ezicl=ezicl+ezicl*rmcrcl den = alb - ezicl pden = palb - ezicl if (den.ne.0.) then tcll=10.0**(tlogcl) tcll=tcll+tcll*rmcrcl ref = 1. / (tcll/den + sbcld) else ref = 0. endif if (pden.ne.0.) then tcll=10.0**(tlogcl) tcll=tcll+tcll*rmcrcl pref = 1. / (tcll/pden + sbcld) else ref = 0. endif clfrac = 1.0 pcfrac = 1.0 clref = ref pcref = pref c if(frstrf) dnd331 = (-100.0*alog10(radcl(2)/radcl(1)))/ 1 (omeghi-omeglo) c c -- calculate reflectivity using cloud fraction and nominal refls c else ref = grref + clfrac*(clref-grref) pref = pgrref + pcfrac*(pclref-pgrref) c if(frstrf) dnd331 = (-100.0*alog10( 1 ((1-clfrac)*radgr(2)+clfrac*(radcl(2)-radgr(2)))/ 2 ((1-clfrac)*radgr(1)+clfrac*(radcl(1)-radgr(1)))))/ 3 (omeghi-omeglo) c endif c dndrm = (pxnvalm-xnvalm(irefl))/(pref-ref) c if (lprint(11)) write (6,1300) pcfrac,pref,clfrac,ref c frstrf = .false. c return c c -- format statements c 1000 format (/'Subroutine reflec'/'Input: ','xlat = ',f8.3, 1 ' isnow = ',i8,' iofset = ',i8/ 2 ' grref = ',f8.4,' clref = ',f8.4,' irefl = ',i8/ 3 ' pcloud = ',f8.4,' pteran = ',f8.4,' alb = ',f8.4,/, 4 ' ilat = ',i3,' sza = ',f6.2,/, 5 ' estozn = ',f7.1,' fteran = ',f6.3,/, 7 ' frstrf = ',l1,' ramcor = ',6(/8f8.3)) 1200 format ('Final: qgcalc = ',f8.6,' qccalc = ',f8.6, 1 ' den = ',f8.4) 1300 format ('Output: pcfrac = ',f8.4, 1 ' pref = ',f8.5,/, 2 9x,'clfrac = ',f8.4,' ref = ',f8.5) c end