#include subroutine radcsw(pint ,h2ommr ,o3mmr ,aermmr ,rh , $ cld ,clwp ,rel ,rei ,fice , $ eccf ,coszrs ,asdir ,asdif ,aldir , $ aldif ,solin ,qrs ,fsns ,fsnt , $ fsds ,fsnsc ,fsntc ,sols ,soll , $ solsd ,solld ,pardif ,pardir, taucl, tauci) !----------------------------------------------------------------------- ! ! Solar radiation code ! ! Basic method is Delta-Eddington as described in: ! ! Briegleb, Bruce P., 1992: Delta-Eddington ! Approximation for Solar Radiation in the NCAR Community Climate Model, ! Journal of Geophysical Research, Vol 97, D7, pp7603-7612). ! ! Two changes to the basic method described above are: (1) the distinction ! between liquid and ice particle clouds, and (2) the addition of an ! aerosol with sulfate radiative properties. ! ! Divides solar spectrum into 18 intervals from 0.2-5.0 micro-meters. ! solar flux fractions specified for each interval. allows for ! seasonally and diurnally varying solar input. Includes molecular, ! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud, ! and surface absorption. Computes delta-eddington reflections and ! transmissions assuming homogeneously mixed layers. Adds the layers ! assuming scattering between layers to be isotropic, and distinguishes ! direct solar beam from scattered radiation. ! ! Longitude loops are broken into 1 or 2 sections, so that only daylight ! (i.e. coszrs > 0) computations are done. ! ! Note that an extra layer above the model top layer is added. ! ! cgs units are used. ! ! Special diagnostic calculation of the clear sky surface and total column ! absorbed flux is also done for cloud forcing diagnostics. ! ! !---------------------------Code history-------------------------------- ! ! Modified March 1995 to add aerosols ! Original version: B. Briegleb ! Standardized: J. Rosinski, June 1992 ! Reviewed: J. Kiehl, B. Briegleb, August 1992 ! Reviewed: J. Kiehl, April 1996 ! Reviewed: B. Briegleb, May 1996 ! !----------------------------------------------------------------------- ! ! $Id: radcsw.F,v 1.2 2001/09/06 17:19:59 sjlin Exp $ ! $Author: sjlin $ ! #include !------------------------------Parameters------------------------------- #include !----------------------------------------------------------------------- real scon ! Solar constant integer nspint ! Num of spctrl intervals across solar spectrum parameter (scon = 1.367e6) parameter ( nspint = 18 ) !------------------------------Commons---------------------------------- #include ! ! Input arguments ! real pint(plond,plevp) ! Interface pressure real h2ommr(plond,plev) ! Specific humidity (h2o mass mix ratio) real o3mmr(plond,plev) ! Ozone mass mixing ratio real aermmr(plond,plev) ! Aerosol mass mixing ratio real rh(plond,plev) ! Relative humidity (fraction) ! real cld(plond,plevp) ! Fractional cloud cover real clwp(plond,plev) ! Layer liquid water path real rel(plond,plev) ! Liquid effective drop size (microns) real rei(plond,plev) ! Ice effective drop size (microns) real fice(plond,plev) ! Fractional ice content within cloud ! real eccf ! Eccentricity factor (1./earth-sun dist ** 2) real coszrs(plond) ! Cosine solar zenith angle real asdir(plond) ! 0.2-0.7 micro-meter srfc alb to direct rad real aldir(plond) ! 0.7-5.0 micro-meter srfc alb to direct rad real asdif(plond) ! 0.2-0.7 micro-meter srfc alb to diffuse rad real aldif(plond) ! 0.7-5.0 micro-meter srfc alb to diffuse rad ! ! Output arguments ! real solin(plond) ! Incident solar flux real qrs(plond,plev) ! Solar heating rate real fsns(plond) ! Surface absorbed solar flux real fsnt(plond) ! Total column absorbed solar flux real fsds(plond) ! Flux Shortwave Downwelling Surface ! real fsnsc(plond) ! Clear sky surface absorbed solar flux real fsntc(plond) ! Clear sky total column absorbed solar flx real sols(plond) ! Direct solar rad incident on surface (< 0.7) real soll(plond) ! Direct solar rad incident on surface (>= 0.7) real solsd(plond) ! Diffuse solar rad incident on surface (< 0.7) real solld(plond) ! Diffuse solar rad incident on surface (>= 0.7) ! JDC ADD real pardif(plond) ! Diffuse photo synthetically active radiation (0.35-0.70 um) real pardir(plond) ! Direct photo synthetically active radiation (0.35-0.70 um) ! !------------------------------Externals-------------------------------- ! integer isrchfgt ! Search for first array element > 0 integer isrchfle ! Search for first array element < 0 ! !---------------------------Local variables----------------------------- ! integer ns ! Spectral loop index integer i ! Longitude loop index integer k ! Level loop index integer n ! Loop index for daylight integer nloop ! Number of daylight loops integer is(2) ! Daytime start indices integer ie(2) ! Daytime end indices integer indxsl ! Index for cloud particle properties ! ! A. Slingo's data for cloud particle radiative properties (from 'A GCM ! Parameterization for the Shortwave Properties of Water Clouds' JAS ! vol. 46 may 1989 pp 1419-1427) ! real abarl(4) ! A coefficient for extinction optical depth real bbarl(4) ! B coefficient for extinction optical depth real cbarl(4) ! C coefficient for single particle scat albedo real dbarl(4) ! D coefficient for single particle scat albedo real ebarl(4) ! E coefficient for asymmetry parameter real fbarl(4) ! F coefficient for asymmetry parameter save abarl, bbarl, cbarl, dbarl, ebarl, fbarl data abarl/ 2.817e-02, 2.682e-02,2.264e-02,1.281e-02/ data bbarl/ 1.305 , 1.346 ,1.454 ,1.641 / data cbarl/-5.62e-08 ,-6.94e-06 ,4.64e-04 ,0.201 / data dbarl/ 1.63e-07 , 2.35e-05 ,1.24e-03 ,7.56e-03 / data ebarl/ 0.829 , 0.794 ,0.754 ,0.826 / data fbarl/ 2.482e-03, 4.226e-03,6.560e-03,4.353e-03/ real abarli ! A coefficient for current spectral interval real bbarli ! B coefficient for current spectral interval real cbarli ! C coefficient for current spectral interval real dbarli ! D coefficient for current spectral interval real ebarli ! E coefficient for current spectral interval real fbarli ! F coefficient for current spectral interval ! ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor ! greater than 20 micro-meters ! ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836) ! real abari(4) ! a coefficient for extinction optical depth real bbari(4) ! b coefficient for extinction optical depth real cbari(4) ! c coefficient for single particle scat albedo real dbari(4) ! d coefficient for single particle scat albedo real ebari(4) ! e coefficient for asymmetry parameter real fbari(4) ! f coefficient for asymmetry parameter save abari, bbari, cbari, dbari, ebari, fbari data abari/ 3.448e-03, 3.448e-03,3.448e-03,3.448e-03/ data bbari/ 2.431 , 2.431 ,2.431 ,2.431 / data cbari/ 1.00e-05 , 1.10e-04 ,1.861e-02,.46658 / data dbari/ 0.0 , 1.405e-05,8.328e-04,2.05e-05 / data ebari/ 0.7661 , 0.7730 ,0.794 ,0.9595 / data fbari/ 5.851e-04, 5.665e-04,7.267e-04,1.076e-04/ real abarii ! A coefficient for current spectral interval real bbarii ! B coefficient for current spectral interval real cbarii ! C coefficient for current spectral interval real dbarii ! D coefficient for current spectral interval real ebarii ! E coefficient for current spectral interval real fbarii ! F coefficient for current spectral interval ! real delta ! Pressure (atmospheres) for stratos. h2o limit real o2mmr ! O2 mass mixing ratio: save delta, o2mmr data delta / 1.70e-3 / data o2mmr / .23143 / ! ! CO2 info: ! real mmwair ! Mean molecular weight of air real mmwco2 ! Mean molecular weight of co2 real co2mmr ! Co2 mass mixing ratio save mmwair, mmwco2 data mmwair / 28.9644 / data mmwco2 / 44.0000 / real albdir(plond) ! Current spc intrvl srf alb to direct rad real albdif(plond) ! Current spc intrvl srf alb to diffuse rad ! ! Next series depends on spectral interval ! real frcsol(nspint) ! Fraction of solar flux in each spectral interval real wavmin(nspint) ! Min wavelength (micro-meters) of interval real wavmax(nspint) ! Max wavelength (micro-meters) of interval real raytau(nspint) ! Rayleigh scattering optical depth real abh2o(nspint) ! Absorption coefficiant for h2o (cm2/g) real abo3 (nspint) ! Absorption coefficiant for o3 (cm2/g) real abco2(nspint) ! Absorption coefficiant for co2 (cm2/g) real abo2 (nspint) ! Absorption coefficiant for o2 (cm2/g) real ph2o(nspint) ! Weight of h2o in spectral interval real pco2(nspint) ! Weight of co2 in spectral interval real po2 (nspint) ! Weight of o2 in spectral interval save frcsol ,wavmin ,wavmax ,raytau ,abh2o ,abo3 , $ abco2 ,abo2 ,ph2o ,pco2 ,po2 data frcsol / .001488, .001389, .001290, .001686, .002877, $ .003869, .026336, .426131, .526861, .526861, $ .526861, .526861, .526861, .526861, .526861, $ .006239, .001834, .001834/ data wavmin / .200, .245, .265, .275, .285, $ .295, .305, .350, .700, .701, $ .701, .701, .701, .702, .702, $ 2.630, 4.160, 4.160/ data wavmax / .245, .265, .275, .285, .295, $ .305, .350, .700, 5.000, 5.000, $ 5.000, 5.000, 5.000, 5.000, 5.000, $ 2.860, 4.550, 4.550/ data raytau / 4.020, 2.180, 1.700, 1.450, 1.250, $ 1.085, 0.730, 0.135, 0.020, .0001, $ .0001, .0001, .0001, .0001, .0001, $ .0001, .0001, .0001/ ! ! Absorption coefficients ! data abh2o / .000, .000, .000, .000, .000, $ .000, .000, .000, .002, .035, $ .377, 1.950, 9.400, 44.600, 190.000, $ .000, .000, .000/ data abo3 / $ 5.370e+04, 13.080e+04, 9.292e+04, 4.530e+04, 1.616e+04, $ 4.441e+03, 1.775e+02, 2.101e+01, .000, .000, $ .000 , .000 , .000 , .000, .000, $ .000 , .000 , .000 / data abco2 / .000, .000, .000, .000, .000, $ .000, .000, .000, .000, .000, $ .000, .000, .000, .000, .000, $ .094, .196, 1.963/ data abo2 / .000, .000, .000, .000, .000, $ .000, .000,1.11e-05,6.69e-05, .000, $ .000, .000, .000, .000, .000, $ .000, .000, .000/ ! ! Spectral interval weights ! data ph2o / .000, .000, .000, .000, .000, $ .000, .000, .000, .505, .210, $ .120, .070, .048, .029, .018, $ .000, .000, .000/ data pco2 / .000, .000, .000, .000, .000, $ .000, .000, .000, .000, .000, $ .000, .000, .000, .000, .000, $ 1.000, .640, .360/ data po2 / .000, .000, .000, .000, .000, $ .000, .000, 1.000, 1.000, .000, $ .000, .000, .000, .000, .000, $ .000, .000, .000/ ! ! Diagnostic and accumulation arrays; note that sfltot, fswup, and ! fswdn are not used in the computation,but are retained for future use. ! real solflx(plond) ! Solar flux in current interval real sfltot(plond) ! Spectrally summed total solar flux real totfld(plond,0:plev) ! Spectrally summed flux divergence real fswup(plond,0:plevp) ! Spectrally summed up flux real fswdn(plond,0:plevp) ! Spectrally summed down flux ! ! Cloud radiative property arrays ! real taucl(plond,plev) ! output cloud water optical depth for visile (0.25-0.69um) ! (no level for above model top) real tauci(plond,plev) ! output cloud ice optical depth for visible ().25-0.69um) ! (no level for above model top) real tauxcl(plond,0:plev) ! water cloud extinction optical depth real tauxci(plond,0:plev) ! ice cloud extinction optical depth real wcl(plond,0:plev) ! liquid cloud single scattering albedo real gcl(plond,0:plev) ! liquid cloud asymmetry parameter real fcl(plond,0:plev) ! liquid cloud forward scattered fraction real wci(plond,0:plev) ! ice cloud single scattering albedo real gci(plond,0:plev) ! ice cloud asymmetry parameter real fci(plond,0:plev) ! ice cloud forward scattered fraction ! ! Aerosol radiative property arrays ! real tauxar(plond,0:plev) ! aerosol extinction optical depth real wa(plond,0:plev) ! aerosol single scattering albedo real ga(plond,0:plev) ! aerosol assymetry parameter real fa(plond,0:plev) ! aerosol forward scattered fraction real tauaer(plond) ! total column aerosol extinction real waer(plond) ! aerosol single scattering albedo real gaer(plond) ! aerosol asymmetry parameter real faer(plond) ! aerosol forward scattering fraction ! ! Sulphate aerosol properties taken from: ! ! Kiehl, J.T., B.P.Briegleb, 1993. The Relative Roles of Sulfate Aerosols ! and Greenhouse Gases in Climate Forcing. Science, Vol. 260, pp. 311-314. ! real ksa(nspint) ! aerosol spectral mass absorption coeff(m2/g) real wsa(nspint) ! aerosol spectral single scattering albedo real gsa(nspint) ! aerosol spectral asymmetry parameter ! data ksa /11.1163, 10.5472, 10.2468, 10.0392, 9.8292, $ 9.6199, 9.0407, 5.3012, 1.9169, 0.3780, $ 0.3780, 0.3780, 0.3780, 0.5704, 0.5704, $ 0.5704, 0.5704, 0.5704 / data wsa / .999999, .999999, .999999, .999999, .999999, $ .999999, .999999, .999999, .999991, .989772, $ .989772, .989772, .989772, .847061, .847061, $ .847061, .847061, .847061 / data gsa / .719161, .719012, .718453, .717820, .716997, $ .715974, .712743, .694889, .618115, .485286, $ .485286, .485286, .485286, .295557, .295557, $ .295557, .295557, .295557 / ! ! Other variables and arrays needed for aerosol: ! real rhfac ! multiplication factor for kaer real rhpc ! level relative humidity in % real a0 ! constant in rh mult factor real a1 ! constant in rh mult factor real a2 ! constant in rh mult factor real a3 ! constant in rh mult factor save a0,a1,a2,a3 data a0 / -9.2906106183 / data a1 / 0.52570211505 / data a2 / -0.0089285760691 / data a3 / 5.0877212432e-05/ ! ! Various arrays and other constants: ! real pflx(plond,0:plevp) ! Interface press, including extra layer real zenfac(plond) ! Square root of cos solar zenith angle real sqrco2 ! Square root of the co2 mass mixg ratio real tmp1 ! Temporary constant array real tmp2 ! Temporary constant array real pdel ! Pressure difference across layer real path ! Mass path of layer real ptop ! Lower interface pressure of extra layer real ptho2 ! Used to compute mass path of o2 real ptho3 ! Used to compute mass path of o3 real pthco2 ! Used to compute mass path of co2 real pthh2o ! Used to compute mass path of h2o real h2ostr ! Inverse square root h2o mass mixing ratio real wavmid ! Spectral interval middle wavelength real trayoslp ! Rayleigh optical depth/standard pressure real tmp1l ! Temporary constant array real tmp2l ! Temporary constant array real tmp3l ! Temporary constant array real tmp1i ! Temporary constant array real tmp2i ! Temporary constant array real tmp3i ! Temporary constant array real rdenom ! Multiple scattering term real psf ! Frac of solar flux in spect interval real gocp ! Gravity/cp ! ! Layer absorber amounts; note that 0 refers to the extra layer added ! above the top model layer ! real uh2o(plond,0:plev) ! Layer absorber amount of h2o real uo3(plond,0:plev) ! Layer absorber amount of o3 real uco2(plond,0:plev) ! Layer absorber amount of co2 real uo2(plond,0:plev) ! Layer absorber amount of o2 real uaer(plond,0:plev) ! Layer aerosol amount ! ! Total column absorber amounts: ! real uth2o(plond) ! Total column absorber amount of h2o real uto3(plond) ! Total column absorber amount of o3 real utco2(plond) ! Total column absorber amount of co2 real uto2(plond) ! Total column absorber amount of o2 real utaer(plond) ! Total column aerosol ! ! These arrays are defined for plev model layers; 0 refers to the extra ! layer on top: ! real rdir(plond,0:plev) ! Layer reflectivity to direct rad real rdif(plond,0:plev) ! Layer reflectivity to diffuse rad real tdir(plond,0:plev) ! Layer transmission to direct rad real tdif(plond,0:plev) ! Layer transmission to diffuse rad real explay(plond,0:plev) ! Solar beam exp transmission for layer real flxdiv(plond,0:plev) ! Flux divergence for layer ! ! These arrays are defined at model interfaces; 0 is the top of the ! extra layer above the model top; plevp is the earth surface: ! real rupdir(plond,0:plevp) ! Ref to dir rad for layers below real rupdif(plond,0:plevp) ! Ref to dif rad for layers below real rdndif(plond,0:plevp) ! Ref to dif rad for layers above real exptdn(plond,0:plevp) ! Solar beam exp down transm from top real tottrn(plond,0:plevp) ! Total transmission for layers above real fluxup(plond,0:plevp) ! Up flux at model interface real fluxdn(plond,0:plevp) ! Down flux at model interface ! !----------------------------------------------------------------------- ! ! Initialize output fields: ! do i=1, plon fsds(i) = 0.0 fsnt(i) = 0.0 fsns(i) = 0.0 solin(i) = 0.0 fsnsc(i) = 0.0 fsntc(i) = 0.0 sols(i) = 0.0 soll(i) = 0.0 solsd(i) = 0.0 solld(i) = 0.0 pardif(i) = 0.0 pardir(i) = 0.0 end do do k=1, plev do i=1, plon qrs(i,k) = 0.0 end do end do ! JDC Added ! ! compute taucl and tauci , even at night ! do k=1,plev do i=1,plon tmp1l = abarl(1) + bbarl(1)/rel(i,k) taucl(i,k) = clwp(i,k)*tmp1l*(1.-fice(i,k)) $ * cld(i,k)*sqrt(cld(i,k)) taucl(i,k)=max(0.0,taucl(i,k)) tmp1i = abari(1) + bbari(1)/rei(i,k) tauci(i,k) = clwp(i,k)*tmp1i*fice(i,k) $ * cld(i,k)*sqrt(cld(i,k)) tauci(i,k)=max(0.0,tauci(i,k)) end do end do ! ! Compute starting, ending daytime loop indices: ! nloop = 0 is(1) = isrchfgt(plon,coszrs,1,0.0) ! ! If night everywhere, return: ! if(is(1).gt.plon) return ie(1) = isrchfle(plon-is(1),coszrs(is(1)+1),1,0.0) + is(1) - 1 nloop = 1 ! ! Possibly 2 daytime loops needed: ! if (ie(1).ne.plon) then is(2) = isrchfgt(plon-ie(1),coszrs(ie(1)+1),1,0.0) + ie(1) if(is(2).le.plon) then nloop = 2 ie(2) = plon end if end if ! ! Define solar incident radiation and interface pressures: ! do n=1,nloop do i=is(n),ie(n) solin(i) = scon*eccf*coszrs(i) pflx(i,0) = 0. end do end do do k=1,plevp do n=1,nloop do i=is(n),ie(n) pflx(i,k) = pint(i,k) end do end do end do ! ! Compute optical paths: ! tmp1 = 0.5/(gravit*sslp) co2mmr = co2vmr*(mmwco2/mmwair) sqrco2 = sqrt(co2mmr) do n=1,nloop do i=is(n),ie(n) ptop = pflx(i,1) ptho2 = o2mmr * ptop / gravit ptho3 = o3mmr(i,1) * ptop / gravit pthco2 = sqrco2 * (ptop / gravit) h2ostr = sqrt( 1. / h2ommr(i,1) ) zenfac(i) = sqrt(coszrs(i)) pthh2o = ptop**2*tmp1 + (ptop*rga)*(h2ostr*zenfac(i)*delta) uh2o(i,0) = h2ommr(i,1)*pthh2o uco2(i,0) = zenfac(i)*pthco2 uo2 (i,0) = zenfac(i)*ptho2 uo3 (i,0) = ptho3 uaer(i,0) = 0.0 end do end do ! tmp2 = delta/gravit do k=1,plev do n=1,nloop do i=is(n),ie(n) pdel = pflx(i,k+1) - pflx(i,k) path = pdel / gravit ptho2 = o2mmr * path ptho3 = o3mmr(i,k) * path pthco2 = sqrco2 * path h2ostr = sqrt(1.0/h2ommr(i,k)) pthh2o = (pflx(i,k+1)**2 - pflx(i,k)**2)*tmp1 + $ pdel*h2ostr*zenfac(i)*tmp2 uh2o(i,k) = h2ommr(i,k)*pthh2o uco2(i,k) = zenfac(i)*pthco2 uo2 (i,k) = zenfac(i)*ptho2 uo3 (i,k) = ptho3 ! ! Adjust aerosol amount by relative humidity factor: ! if( rh(i,k) .gt. .90 ) then rhfac = 2.8 else if (rh(i,k) .lt. .60 ) then rhfac = 1.0 else rhpc = 100. * rh(i,k) rhfac = (a0 + a1*rhpc + a2*rhpc**2 + a3*rhpc**3) endif uaer(i,k) = aermmr(i,k)*rhfac*path end do end do end do ! ! Compute column absorber amounts for the clear sky computation: ! do n=1,nloop do i=is(n),ie(n) uth2o(i) = 0.0 uto3(i) = 0.0 utco2(i) = 0.0 uto2(i) = 0.0 utaer(i) = 0.0 end do end do do k=1,plev do n=1,nloop do i=is(n),ie(n) uth2o(i) = uth2o(i) + uh2o(i,k) uto3(i) = uto3(i) + uo3(i,k) utco2(i) = utco2(i) + uco2(i,k) uto2(i) = uto2(i) + uo2(i,k) utaer(i) = utaer(i) + uaer(i,k) end do end do end do ! ! Initialize spectrally integrated totals: ! do k=0,plev do i=1,plon totfld(i,k) = 0.0 fswup (i,k) = 0.0 fswdn (i,k) = 0.0 end do end do do i=1,plon sfltot(i) = 0.0 fswup (i,plevp) = 0.0 fswdn (i,plevp) = 0.0 end do ! ! Set cloud properties for top (0) layer; so long as tauxcl is zero, ! there is no cloud above top of model; the other cloud properties ! are arbitrary: ! do n=1,nloop do i=is(n),ie(n) tauxcl(i,0) = 0. wcl(i,0) = 0.999999 gcl(i,0) = 0.85 fcl(i,0) = 0.725 tauxci(i,0) = 0. wci(i,0) = 0.999999 gci(i,0) = 0.85 fci(i,0) = 0.725 ! ! Aerosol ! tauxar(i,0) = 0. wa(i,0) = 0.925 ga(i,0) = 0.850 fa(i,0) = 0.7225 end do end do ! ! Begin spectral loop ! do 100 ns=1,nspint ! ! Set index for cloud particle properties based on the wavelength, ! according to A. Slingo (1989) equations 1-3: ! Use index 1 (0.25 to 0.69 micrometers) for visible ! Use index 2 (0.69 - 1.19 micrometers) for near-infrared ! Use index 3 (1.19 to 2.38 micrometers) for near-infrared ! Use index 4 (2.38 to 4.00 micrometers) for near-infrared ! ! Note that the minimum wavelength is encoded (with .001, .002, .003) ! in order to specify the index appropriate for the near-infrared ! cloud absorption properties ! if(wavmax(ns) .le. 0.7) then indxsl = 1 else if(wavmin(ns) .eq. 0.700) then indxsl = 2 else if(wavmin(ns) .eq. 0.701) then indxsl = 3 else if(wavmin(ns) .eq. 0.702 .or. wavmin(ns) .gt. 2.38) then indxsl = 4 end if ! ! Set cloud extinction optical depth, single scatter albedo, ! asymmetry parameter, and forward scattered fraction: ! abarli = abarl(indxsl) bbarli = bbarl(indxsl) cbarli = cbarl(indxsl) dbarli = dbarl(indxsl) ebarli = ebarl(indxsl) fbarli = fbarl(indxsl) abarii = abari(indxsl) bbarii = bbari(indxsl) cbarii = cbari(indxsl) dbarii = dbari(indxsl) ebarii = ebari(indxsl) fbarii = fbari(indxsl) do k=1,plev do n=1,nloop do i=is(n),ie(n) ! ! liquid ! tmp1l = abarli + bbarli/rel(i,k) tmp2l = 1. - cbarli - dbarli*rel(i,k) tmp3l = fbarli*rel(i,k) ! ! ice ! tmp1i = abarii + bbarii/rei(i,k) tmp2i = 1. - cbarii - dbarii*rei(i,k) tmp3i = fbarii*rei(i,k) ! ! Cloud fraction incorporated into cloud extinction optical depth ! tauxcl(i,k) = clwp(i,k)*tmp1l*(1.-fice(i,k)) $ *cld(i,k)*sqrt(cld(i,k)) tauxci(i,k) = clwp(i,k)*tmp1i*fice(i,k) $ *cld(i,k)*sqrt(cld(i,k)) ! ! Do not let single scatter albedo be 1; delta-eddington solution ! for non-conservative case: ! wcl(i,k) = min(tmp2l,.999999) gcl(i,k) = ebarli + tmp3l fcl(i,k) = gcl(i,k)*gcl(i,k) wci(i,k) = min(tmp2i,.999999) gci(i,k) = ebarii + tmp3i fci(i,k) = gci(i,k)*gci(i,k) ! ! Set aerosol properties ! Conversion factor to adjust aerosol extinction (m2/g) ! tauxar(i,k) = 1.e4 * ksa(ns) * uaer(i,k) wa(i,k) = wsa(ns) ga(i,k) = gsa(ns) fa(i,k) = gsa(ns)*gsa(ns) waer(i) = wa(i,k) gaer(i) = ga(i,k) faer(i) = fa(i,k) end do end do end do ! ! Set reflectivities for surface based on mid-point wavelength ! wavmid = 0.5*(wavmin(ns) + wavmax(ns)) ! ! Wavelength less than 0.7 micro-meter ! if (wavmid .lt. 0.7 ) then do n=1,nloop do i=is(n),ie(n) albdir(i) = asdir(i) albdif(i) = asdif(i) end do end do ! ! Wavelength greater than 0.7 micro-meter ! else do n=1,nloop do i=is(n),ie(n) albdir(i) = aldir(i) albdif(i) = aldif(i) end do end do end if trayoslp = raytau(ns)/sslp ! ! Layer input properties now completely specified; compute the ! delta-Eddington solution reflectivities and transmissivities ! for each layer, starting from the top and working downwards: ! call t_startf('radded') call radded(coszrs ,trayoslp,pflx ,abh2o(ns),abo3(ns), $ abco2(ns),abo2(ns),uh2o ,uo3 ,uco2 , $ uo2 ,tauxcl ,wcl ,gcl ,fcl , $ tauxci ,wci ,gci ,fci ,tauxar , $ wa ,ga ,fa ,nloop ,is , $ ie ,rdir ,rdif ,tdir ,tdif , $ explay ,exptdn ,rdndif ,tottrn ) call t_stopf('radded') ! ! Compute reflectivity to direct and diffuse radiation for layers below ! by adding succesive layers starting from the surface and working ! upwards: ! do n=1,nloop do i=is(n),ie(n) rupdir(i,plevp) = albdir(i) rupdif(i,plevp) = albdif(i) end do end do do k=plev,0,-1 do n=1,nloop do i=is(n),ie(n) rdenom = 1./( 1. - rdif(i,k)*rupdif(i,k+1)) rupdir(i,k) = rdir(i,k) + tdif(i,k)* $ (rupdir(i,k+1)*explay(i,k) + $ rupdif(i,k+1)*(tdir(i,k)-explay(i,k)))*rdenom rupdif(i,k) = rdif(i,k) + $ rupdif(i,k+1)*tdif(i,k)**2*rdenom end do end do end do ! ! Compute up and down fluxes for each interface, using the added ! atmospheric layer properties at each interface: ! do k=0,plevp do n=1,nloop do i=is(n),ie(n) rdenom = 1./(1. - rdndif(i,k)*rupdif(i,k)) fluxup(i,k) = (exptdn(i,k)*rupdir(i,k) + $ (tottrn(i,k)-exptdn(i,k))*rupdif(i,k))*rdenom fluxdn(i,k)=exptdn(i,k) + (tottrn(i,k) - exptdn(i,k) + $ exptdn(i,k)*rupdir(i,k)*rdndif(i,k))*rdenom end do end do end do ! ! Compute flux divergence in each layer using the interface up and down ! fluxes: ! do k=0,plev do n=1,nloop do i=is(n),ie(n) flxdiv(i,k) = (fluxdn(i,k ) - fluxdn(i,k+1)) + $ (fluxup(i,k+1) - fluxup(i,k )) end do end do end do ! ! Monochromatic computation completed; accumulate in totals; adjust ! fraction within spectral interval to allow for the possibility of ! sub-divisions within a particular interval: ! psf = 1.0 if(ph2o(ns).ne.0.) psf = psf*ph2o(ns) if(pco2(ns).ne.0.) psf = psf*pco2(ns) if(po2 (ns).ne.0.) psf = psf*po2 (ns) do n=1,nloop do i=is(n),ie(n) solflx(i) = solin(i)*frcsol(ns)*psf fsnt(i) = fsnt(i) + solflx(i)*(fluxdn(i,1) - fluxup(i,1)) fsns(i) = fsns(i) + solflx(i)* $ (fluxdn(i,plevp) - fluxup(i,plevp)) sfltot(i) = sfltot(i) + solflx(i) fswup(i,0) = fswup(i,0) + solflx(i)*fluxup(i,0) fswdn(i,0) = fswdn(i,0) + solflx(i)*fluxdn(i,0) ! ! Down spectral fluxes need to be in mks; thus the .001 conversion factors ! if (wavmid .lt. 0.7) then sols(i) = sols(i) + exptdn(i,plevp)*solflx(i)*0.001 solsd(i) = solsd(i) + (fluxdn(i,plevp) - $ exptdn(i,plevp)) * solflx(i)*0.001 if(wavmid .gt. 0.35) then pardir(i) = pardir(i) + exptdn(i,plevp)*solflx(i)*0.001 pardif(i) = pardif(i) + (fluxdn(i,plevp) - $ exptdn(i,plevp)) * solflx(i)*0.001 end if else soll(i) = soll(i) + exptdn(i,plevp)*solflx(i)*0.001 solld(i) = solld(i) + (fluxdn(i,plevp) - $ exptdn(i,plevp)) * solflx(i)*0.001 end if end do end do do k=0,plev do n=1,nloop do i=is(n),ie(n) totfld(i,k) = totfld(i,k) + solflx(i)*flxdiv(i,k) fswup(i,k+1) = fswup(i,k+1) + solflx(i)*fluxup(i,k+1) fswdn(i,k+1) = fswdn(i,k+1) + solflx(i)*fluxdn(i,k+1) end do end do end do ! ! Following code is the diagnostic clear sky computation: ! ! Compute delta-Eddington solution reflectivities and transmissivities ! for the entire column; note, for convenience, we use the same ! reflectivity and transmissivity arrays as for the full calculation ! above, where 0 for layer quantities refers to the entire atmospheric ! column, and where 0 for interface quantities refers to top of atmos- ! phere, while 1 refers to the surface: ! ! ! Compute total column aerosol optical depth: ! do n=1,nloop do i=is(n),ie(n) ! ! Conversion factor to adjust aerosol extinction (m2/g) ! tauaer(i) = 1.e4 * ksa(ns) * utaer(i) end do end do call radclr(coszrs ,trayoslp,pflx ,abh2o(ns),abo3(ns) , $ abco2(ns),abo2(ns),uth2o ,uto3 ,utco2 , $ uto2 ,tauaer ,waer ,gaer ,faer , $ nloop ,is ,ie ,rdir ,rdif , $ tdir ,tdif ,explay ,exptdn ,rdndif , $ tottrn ) ! ! Compute reflectivity to direct and diffuse radiation for entire ! column; 0,1 on layer quantities refers to two effective layers ! overlying surface; 0 on interface quantities refers to top of column; ! 2 on interface quantities refers to the surface: ! do n=1,nloop do i=is(n),ie(n) rupdir(i,2) = albdir(i) rupdif(i,2) = albdif(i) end do end do do k=1,0,-1 do n=1,nloop do i=is(n),ie(n) rdenom = 1./( 1. - rdif(i,k)*rupdif(i,k+1)) rupdir(i,k) = rdir(i,k) + tdif(i,k)* $ (rupdir(i,k+1)*explay(i,k) + $ rupdif(i,k+1)*(tdir(i,k)-explay(i,k)))*rdenom rupdif(i,k) = rdif(i,k) + $ rupdif(i,k+1)*tdif(i,k)**2*rdenom end do end do end do ! ! Compute up and down fluxes for each interface, using the added ! atmospheric layer properties at each interface: ! do k=0,2 do n=1,nloop do i=is(n),ie(n) rdenom = 1./(1. - rdndif(i,k)*rupdif(i,k)) fluxup(i,k) = (exptdn(i,k)*rupdir(i,k) + $ (tottrn(i,k)-exptdn(i,k))*rupdif(i,k))*rdenom fluxdn(i,k)=exptdn(i,k) + (tottrn(i,k) - exptdn(i,k) + $ exptdn(i,k)*rupdir(i,k)*rdndif(i,k))*rdenom end do end do end do do n=1,nloop do i=is(n),ie(n) fsntc(i) = fsntc(i) + solflx(i)*(fluxdn(i,0)-fluxup(i,0)) fsnsc(i) = fsnsc(i) + solflx(i)*(fluxdn(i,2)-fluxup(i,2)) end do end do ! ! End of clear sky calculation ! 100 continue ! End of spectral interval loop ! ! Compute solar heating rate (k/s) ! gocp = gravit/cpair do k=1,plev do n=1,nloop do i=is(n),ie(n) qrs(i,k) = -gocp*totfld(i,k)/(pint(i,k) - pint(i,k+1)) end do end do end do ! ! Set the downwelling flux at the surface ! do i=1,plon fsds(i) = fswdn(i,plevp) end do return end