SUBROUTINE tuv(nz, z, airlev, alb, coszen, $ odcld, omcld, gcld, $ odaer1, odaer2, omaer1, omaer2, gaer, $ rafcld) IMPLICIT NONE **** Inputs: * altitude grid INTEGER nz REAL z(nz) * solar zenith angle REAL coszen * air profile REAL airlev(nz) * cloud and aerosol REAL odcld(nz), omcld(nz), gcld(nz) REAL odaer1(nz), omaer1(nz), gaer(nz) REAL odaer2(nz), omaer2(nz) * surface albedo REAL alb **** Output: Actinic flux cloudy/clear ratio as a funct of altitude REAL rafcld(nz) **** Internal: * internally computed optical depths (to avoid overwriting odcld, odaer) REAL dtrl(nz), dtcld(nz), dtaer(nz) * Rayleigh scattering cross section from WMO 1985 (originally from * Nicolet, M., On the molecular scattering in the terrestrial atmosphere: * An empirical formula for its calculation in the homoshpere, Planet. * Space Sci., 32, 1467-1468, 1984. * computed at the selected wavelenth wc(nm) using the equations: c wmicrn = wc/1.E3 c IF( wmicrn .LE. 0.55) THEN c xx = 3.6772 + 0.389*wmicrn + 0.09426/wmicrn c ELSE c xx = 4. + 0.04 c ENDIF c srayl(iw) = 4.02e-28/(wmicrn)**xx REAL srayl DATA srayl /2.922E-26/ * scale heights for air density REAL hscale DATA hscale /8.05e5/ * actinic flux (relative units) REAL fdir(nz), fdn(nz), fup(nz) REAL actflx(2,nz) * loop indices INTEGER i * ___ * compute Rayleigh optical depths: * include exponential tail integral at top (iz = nz) DO i = 1, nz-1 dtrl(i) = (airlev(i+1)+airlev(i)) * 0.5 * $ (z(i+1)-z(i)) * 1.E5 * srayl ENDDO dtrl(nz) = hscale * airlev(nz) * 1.E5 **** clear sky calculation ** monochromatic radiative transfer: * call rtlink for plane-parallel 2-stream DO i = 1, nz-1 dtcld(i) = 0. dtaer(i) = odaer1(i) ENDDO CALL rtlpp2(nz, alb, coszen, $ dtrl, $ dtcld, omcld, gcld, $ dtaer,omaer1,gaer, $ fdir, fdn, fup) DO i = 1, nz actflx(1,i) = fdir(i) + fdn(i) + fup(i) ENDDO **** cloudy sky calculation: DO i = 1, nz-1 dtcld(i) = odcld(i) dtaer(i) = odaer2(i) ENDDO CALL rtlpp2(nz, alb, coszen, $ dtrl, $ dtcld, omcld, gcld, $ dtaer,omaer2,gaer, $ fdir, fdn, fup) DO i = 1, nz actflx(2,i) = fdir(i) + fdn(i) + fup(i) ENDDO DO i = 1, nz rafcld(i) = actflx(2,i)/actflx(1,i) ENDDO RETURN END *============================================================================== SUBROUTINE rtlpp2(nz, ag, m0, $ dtrl, $ dtcld, omcld, gcld, $ dtaer,omaer,gaer, $ fdir, fdn, fup) *_______________________________________________________________________ IMPLICIT NONE * input: INTEGER nz REAL ag REAL m0 REAL dtrl(nz) REAL dtcld(nz), omcld(nz), gcld(nz) REAL dtaer(nz), omaer(nz), gaer(nz) * output REAL fdir(nz), fdn(nz), fup(nz) * local: REAL dt(nz), om(nz), g(nz) REAL fdiri(nz), fdni(nz), fupi(nz) REAL daaer, dtsct, dtabs, dsaer, dscld, dacld INTEGER i, ii LOGICAL delta DATA delta /.true./ * largest number of the machine: REAL largest PARAMETER(largest=1.E+36) *_______________________________________________________________________ * initialize: DO 5 i = 1, nz fdir(i) = 0. fup(i) = 0. fdn(i) = 0. 5 CONTINUE * ---- * set here any coefficients specific to rt scheme, * ---- DO 10, i = 1, nz - 1 dscld = dtcld(i)*omcld(i) dacld = dtcld(i)*(1.-omcld(i)) dsaer = dtaer(i)*omaer(i) daaer = dtaer(i)*(1.-omaer(i)) dtsct = dtrl(i) + dscld + dsaer dtabs = dacld + daaer dtabs = AMAX1(dtabs,1./largest) dtsct = AMAX1(dtsct,1./largest) * invert z-coordinate: ii = nz - i dt(ii) = dtsct + dtabs om(ii) = dtsct/(dtsct + dtabs) IF(dtsct .EQ. 1./largest) om(ii) = 1./largest g(ii) = (gcld(i)*dscld + gaer(i)*dsaer)/dtsct 10 CONTINUE * call rt routine: CALL pp2str(nz,m0,ag,dt,om,g, $ fdiri, fupi, fdni) * put on upright z-coordinate DO 20, i = 1, nz ii = nz - i + 1 fdir(i) = fdiri(ii) fup(i) = fupi(ii) fdn(i) = fdni(ii) 20 CONTINUE *_______________________________________________________________________ RETURN END *============================================================================== SUBROUTINE pp2str(nz,mu,rsfc,tauu,omu,gu, $ fdr, fup, fdn) *_______________________________________________________________________ * TWO-STREAM EQUATIONS FOR MULTIPLE LAYERS * based on equations from Toon et al., Journal of Geophysical Research * Volume 94, #D13 Nov. 20, 1989 Issue * programmed by: Kathleen G. Mosher * for Sasha Madronich, N.C.A.R. A.C.D. *_______________________________________________________________________ * Now it contains 9 two-stream methods to choose from. * programmed on 05.26.94 by: Irina V.Petropavlovskikh * for Sasha Madronich, N.C.A.R. A.C.D. *_______________________________________________________________________ IMPLICIT NONE ******* * input: ******* INTEGER nz REAL mu, rsfc REAL tauu(nz), omu(nz), gu(nz) ******* * output: ******* REAL fup(nz),fdn(nz),fdr(nz) REAL eup(nz),edn(nz),edr(nz) ******* * local: ******* * internal coefficients and matrix REAL lam(nz),taun(nz),bgam(nz) REAL e1(nz),e2(nz),e3(nz),e4(nz) REAL cup(nz),cdn(nz),cuptn(nz),cdntn(nz) REAL tauc, mu1(nz) INTEGER row REAL a(2*nz),b(2*nz),d(2*nz),e(2*nz),y(2*nz) ******* * other: ******* REAL pifs, fdn0 REAL f, g, tau, om, tempg REAL gam1, gam2, gam3, gam4 * For calculations of Associated Legendre Polynomials for GAMA1,2,3,4 * in delta-function, modified quadrature, hemispheric constant, * Hybrid modified Eddington-delta function metods, p633,Table1. * W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 * W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, * uncomment the following two lines and the appropriate statements further * down. C REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0, C > BETA1, BETAn, amu1, subd REAL expon, expon0, expon1, divisr, temp, up, dn REAL ssfc REAL taug INTEGER nlayer, mrows, lev INTEGER i, j * Some additional program constants: REAL eps, precis PARAMETER (eps = 1.E-3, precis = 1.E-7) *_______________________________________________________________________ * MU = cosine of solar zenith angle * RSFC = surface albedo * TAUU = unscaled optical depth of each layer * OMU = unscaled single scattering albedo * GU = unscaled asymmetry factor * KLEV = max dimension of number of layers in atmosphere * NLAYER = number of layers in the atmosphere * NLEVEL = nlayer + 1 = number of levels * initial conditions: pi*solar flux = 1; diffuse incidence = 0 pifs = 1. fdn0 = 0. nlayer = nz - 1 ************** compute coefficients for each layer: * GAM1 - GAM4 = 2-stream coefficients, different for different approximations * EXPON0 = calculation of e when TAU is zero * EXPON1 = calculation of e when TAU is TAUN * CUP and CDN = calculation when TAU is zero * CUPTN and CDNTN = calc. when TAU is TAUN * DIVISR = prevents division by zero tauc = 0. DO 10, i = 1, nlayer g = gu(i) tau = tauu(i) om = omu(i) * stay away from 1 by precision. For g, also stay away from -1 tempg = AMIN1(abs(g),1. - precis) g = SIGN(tempg,g) om = AMIN1(om,1-precis) * delta-scaling. Have to be done for delta-Eddington approximation, * delta discrete ordinate, Practical Improved Flux Method, delta function, * and Hybrid modified Eddington-delta function methods approximations f = g*g g = (g - f)/(1 - f) taun(i) = (1 - om*f)*tau om = (1 - f)*om/(1 - om*f) *** the following gamma equations are from pg 16,289, Table 1 * Eddington approximation(Joseph et al., 1976, JAS, 33, 2452): gam1 = (7. - om*(4. + 3.*g))/4. gam2 = -(1. - om*(4. - 3.*g))/4. gam3 = (2. - 3.*g*mu)/4. gam4 = 1. - gam3 mu1(i) = 0.5 * quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475): c gam1 = 1.7320508*(2. - om*(1. + g))/.2 c gam2 = 1.7320508*om*(1. - g)/2. c gam3 = (1. - 1.7320508*g*mu)/2 c gam4 = 1. - gam3 * hemispheric mean (Toon et al., 1089, JGR, 94, 16287): c gam1 = 2. - om*(1. + g) c gam2 = om*(1. - g) c gam3 = (2. - g*mu)/4. c gam4 = 1. - gam3 * PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166): c GAM1 = 0.25*(8. - OM*(5. + 3.*G)) c GAM2 = 0.75*OM*(1.-G) c GAM3 = 0.25*(2.-3.*G*MU) c GAM4 = 1. - GAM3 * delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26): c GAM1 = 0.5*1.7320508*(2. - OM*(1. + G)) c GAM2 = 0.5*1.7320508*OM*(1.-G) c GAM3 = 0.5*(1.-1.7320508*G*MU) c GAM4 = 1. - GAM3 * Calculations of Associated Legendre Polynomials for GAMA1,2,3,4 * in delta-function, modified quadrature, hemispheric constant, * Hybrid modified Eddington-delta function metods, p633,Table1. * W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 * W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440 c YLM0 = 2. c YLM2 = -3.*G*MU c YLM4 = 0.875*G**3*MU*(5.*MU**2-3.) c YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4) c YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4 c *+429.*MU**6) c YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4 c *-25740.*MU**6+12155.*MU**8) c YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4 c *+218790.*MU**6-230945.*MU**8+88179.*MU**10) c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 c YLMS=0.25*YLMS c BETA0 = YLMS c c amu1=1./1.7320508 c YLM0 = 2. c YLM2 = -3.*G*amu1 c YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.) c YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4) c YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4 c *+429.*amu1**6) c YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4 c *-25740.*amu1**6+12155.*amu1**8) c YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4 c *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10) c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 c YLMS=0.25*YLMS c BETA1 = YLMS c c BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5 c *-0.045776*G**7) * Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630): c subd=4.*(1.-G*G*(1.-MU)) c GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd c GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd c GAM3 = BETA0 c GAM4 = 1. - GAM3 ***** * delta function (Meador, and Weaver, 1980, JAS, 37, 630): c GAM1 = (1. - OM*(1. - beta0))/MU c GAM2 = OM*BETA0/MU c GAM3 = BETA0 c GAM4 = 1. - GAM3 ***** * modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630): c GAM1 = 1.7320508*(1. - OM*(1. - beta1)) c GAM2 = 1.7320508*OM*beta1 c GAM3 = BETA0 c GAM4 = 1. - GAM3 * hemispheric constant (Toon et al., 1989, JGR, 94, 16287): c GAM1 = 2.*(1. - OM*(1. - betan)) c GAM2 = 2.*OM*BETAn c GAM3 = BETA0 c GAM4 = 1. - GAM3 ***** * save mu1 for use in converting irradiance to actinic flux c mu1(i) = (1-om)/(gam1 - gam2) * lambda = pg 16,290 equation 21 * big gamma = pg 16,290 equation 22 lam(i) = sqrt(gam1*gam1 - gam2*gam2) IF( gam2 .NE. 0.) THEN bgam(i) = (gam1 - lam(i))/gam2 ELSE bgam(i) = 0. ENDIF expon = EXP(-lam(i)*taun(i)) * e1 - e4 = pg 16,292 equation 44 e1(i) = 1. + bgam(i)*expon e2(i) = 1. - bgam(i)*expon e3(i) = bgam(i) + expon e4(i) = bgam(i) - expon * the following sets up for the C equations 23, and 24 * found on page 16,290 * prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3 * which is approx equiv to shifting MU by 0.5*EPS* (MU)**3 expon0 = EXP(-(tauc )/mu) expon1 = EXP(-(tauc + taun(i))/mu) divisr = lam(i)*lam(i) - 1./(mu*mu) temp = AMAX1(eps,abs(divisr)) divisr = SIGN(temp,divisr) up = om*pifs*((gam1 - 1./mu)*gam3 + gam4*gam2)/divisr dn = om*pifs*((gam1 + 1./mu)*gam4 + gam2*gam3)/divisr * cup and cdn are when tau is equal to zero * cuptn and cdntn are when tau is equal to taun cup(i) = up*expon0 cdn(i) = dn*expon0 cuptn(i) = up*expon1 cdntn(i) = dn*expon1 tauc = tauc + taun(i) 10 CONTINUE ***************** set up matrix ****** * ssfc = pg 16,292 equation 37 where pi Fs is one (unity). ssfc = rsfc*mu*EXP(-tauc/mu)*pifs * MROWS = the number of rows in the matrix mrows = 2*nlayer * the following are from pg 16,292 equations 39 - 43. * set up first row of matrix: i = 1 a(1) = 0. b(1) = e1(i) d(1) = -e2(i) e(1) = fdn0 - cdn(i) row=1 * set up odd rows 3 thru (MROWS - 1): i = 0 DO 20, row = 3, mrows - 1, 2 i = i + 1 a(row) = e2(i)*e3(i) - e4(i)*e1(i) b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1) d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1) e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + $ e1(i)*(cdntn(i) - cdn(i + 1)) 20 CONTINUE * set up even rows 2 thru (MROWS - 2): i = 0 DO 30, row = 2, mrows - 2, 2 i = i + 1 a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1) b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1) d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1) e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - $ (cdn(i + 1) - cdntn(i))*e4(i + 1) 30 CONTINUE * set up last row of matrix at MROWS: row = mrows i = nlayer a(row) = e1(i) - rsfc*e3(i) b(row) = e2(i) - rsfc*e4(i) d(row) = 0. e(row) = ssfc - cuptn(i) + rsfc*cdntn(i) * solve tri-diagonal matrix: CALL tridag(a, b, d, e, y, mrows) **** unfold solution of matrix, compute output fluxes: row = 1 lev = 1 j = 1 taug = 0. * the following equations are from pg 16,291 equations 31 & 32 fdr(lev) = 1. edr(lev) = mu edn(lev) = fdn0 eup(lev) = y(row)*e3(j) - y(row + 1)*e4(j) + cup(j) fdn(lev) = edn(lev)/mu1(lev) fup(lev) = eup(lev)/mu1(lev) DO 60, lev = 2, nlayer + 1 taug = taug + taun(j) fdr(lev) = EXP(-taug/mu) edr(lev) = mu*fdr(lev) edn(lev) = y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j) eup(lev) = y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j) fdn(lev) = edn(lev)/mu1(j) fup(lev) = eup(lev)/mu1(j) row = row + 2 j = j + 1 60 CONTINUE *_______________________________________________________________________ RETURN END *============================================================================== SUBROUTINE tridag(a,b,c,r,u,n) *_______________________________________________________________________ * solves tridiagonal system. From Numerical Recipies, p. 40 *_______________________________________________________________________ IMPLICIT NONE * input: INTEGER n REAL a, b, c, r DIMENSION a(n),b(n),c(n),r(n) * output: REAL u DIMENSION u(n) * local: INTEGER j REAL bet, gam DIMENSION gam(n) *_______________________________________________________________________ IF (b(1) .EQ. 0.) STOP 1001 bet = b(1) u(1) = r(1)/bet DO 11, j = 2, n gam(j) = c(j - 1)/bet bet = b(j) - a(j)*gam(j) IF (bet .EQ. 0.) STOP 2002 u(j) = (r(j) - a(j)*u(j - 1))/bet 11 CONTINUE DO 12, j = n - 1, 1, -1 u(j) = u(j) - gam(j + 1)*u(j + 1) 12 CONTINUE *_______________________________________________________________________ RETURN END