!============================================================================= ! This file contains the following subroutines, related to the solution of ! the equation of radiative transfer in multiple homogeneous layers. ! rtlink ! ps2str ! tridag !============================================================================= MODULE rad_trans IMPLICIT NONE private public :: rtlink CONTAINS SUBROUTINE rtlink( & nstr, nlev, nlyr, nwave, & iw, albedo, zen, & dsdh, nid, & dtrl, & dto3, & dto2, & dtso2, & dtno2, & dtcld, omcld, gcld, & dtaer, omaer, gaer, & dtsnw, omsnw, gsnw, & #ifdef SW_DEBUG edir, edn, eup, fdir, fdn, fup, tuv_diags ) #else edir, edn, eup, fdir, fdn, fup ) #endif use params_mod, only : largest, pi !--------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------- INTEGER, intent(in) :: nstr INTEGER, intent(in) :: nlev, nlyr INTEGER, intent(in) :: nwave, iw REAL, intent(in) :: albedo REAL, intent(in) :: zen REAL, intent(in) :: dsdh(0:nlyr,nlyr) INTEGER, intent(in) :: nid(0:nlyr) REAL, intent(in) :: dtrl(nlyr,nwave) REAL, intent(in) :: dto3(nlyr,nwave), dto2(nlyr,nwave) REAL, intent(in) :: dtso2(nlyr,nwave), dtno2(nlyr,nwave) REAL, intent(in) :: dtcld(nlyr,nwave), omcld(nlyr,nwave), gcld(nlyr,nwave) REAL, intent(in) :: dtaer(nlyr,nwave), omaer(nlyr,nwave), gaer(nlyr,nwave) REAL, intent(in) :: dtsnw(nlyr,nwave), omsnw(nlyr,nwave), gsnw(nlyr,nwave) REAL, intent(out) :: edir(nlev), edn(nlev), eup(nlev) REAL, intent(out) :: fdir(nlev), fdn(nlev), fup(nlev) #ifdef SW_DEBUG LOGICAL, intent(in) :: tuv_diags #endif !--------------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------------- REAL, PARAMETER :: dr = pi/180. INTEGER :: k, kk REAL :: dtabs,dtsct,dscld,dsaer,dssnw,dacld,daaer,dasnw REAL :: dt(nlyr), om(nlyr), g(nlyr) !--------------------------------------------------------------------- ! ... specific to ps2str !--------------------------------------------------------------------- LOGICAL, parameter :: delta = .true. REAL ediri(nlev), edni(nlev), eupi(nlev) REAL fdiri(nlev), fdni(nlev), fupi(nlev) !--------------------------------------------------------------------- ! initialize: !--------------------------------------------------------------------- fdir(1:nlev) = 0. fup(1:nlev) = 0. fdn(1:nlev) = 0. edir(1:nlev) = 0. eup(1:nlev) = 0. edn(1:nlev) = 0. DO k = 1, nlyr dscld = dtcld(k,iw)*omcld(k,iw) dacld = dtcld(k,iw)*(1.-omcld(k,iw)) dsaer = dtaer(k,iw)*omaer(k,iw) daaer = dtaer(k,iw)*(1.-omaer(k,iw)) dssnw = dtsnw(k,iw)*omsnw(k,iw) dasnw = dtsnw(k,iw)*(1.-omsnw(k,iw)) dtsct = dtrl(k,iw) + dscld + dsaer + dssnw dtabs = dtso2(k,iw) + dto2(k,iw) + dto3(k,iw) & + dtno2(k,iw) + dacld + daaer + dasnw dtabs = max( dtabs,1./largest ) dtsct = max( dtsct,1./largest ) !--------------------------------------------------------------------- ! from bottom-up to top-down !--------------------------------------------------------------------- kk = nlyr - k + 1 dt(kk) = dtsct + dtabs g(kk) = (gcld(k,iw)*dscld + gsnw(k,iw)*dssnw + gaer(k,iw)*dsaer)/dtsct IF( dtsct /= 1./largest ) then om(kk) = dtsct/(dtsct + dtabs) ELSE om(kk) = 1./largest ENDIF END DO #ifdef SW_DEBUG if( tuv_diags .and. iw == 100 ) then open(unit=33,file='WRF-TUV.dbg1.out') do k = 1,nlyr,5 write(33,'(1p,5g15.7)') dt(k:min(k+4,nlyr)) end do do k = 1,nlyr,5 write(33,'(1p,5g15.7)') g(k:min(k+4,nlyr)) end do do k = 1,nlyr,5 write(33,'(1p,5g15.7)') om(k:min(k+4,nlyr)) end do close(33) endif #endif CALL ps2str( nlyr, nlev, zen, albedo, & dt, om, g, & dsdh, nid, delta, & fdiri, fupi, fdni, ediri, eupi, edni) !--------------------------------------------------------------------- ! from top-down to bottom-up !--------------------------------------------------------------------- fdir(1:nlev) = fdiri(nlev:1:-1) fup(1:nlev) = fupi(nlev:1:-1) fdn(1:nlev) = fdni(nlev:1:-1) edir(1:nlev) = ediri(nlev:1:-1) eup(1:nlev) = eupi(nlev:1:-1) edn(1:nlev) = edni(nlev:1:-1) END SUBROUTINE rtlink SUBROUTINE ps2str( & nlyr, nlev, zen, rsfc, & tauu, omu, gu, & dsdh, nid, delta, & fdr, fup, fdn, edr, eup, edn) !-----------------------------------------------------------------------------* != PURPOSE: =* != Solve two-stream equations for multiple layers. The subroutine is based =* != on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=* != It contains 9 two-stream methods to choose from. A pseudo-spherical =* != correction has also been added. =* !-----------------------------------------------------------------------------* != PARAMETERS: =* != NLEVEL - INTEGER, number of specified altitude levels in the working (I)=* != grid =* != ZEN - REAL, solar zenith angle (degrees) (I)=* != RSFC - REAL, surface albedo at current wavelength (I)=* != TAUU - REAL, unscaled optical depth of each layer (I)=* != OMU - REAL, unscaled single scattering albedo of each layer (I)=* != GU - REAL, unscaled asymmetry parameter of each layer (I)=* != DSDH - REAL, slant path of direct beam through each layer crossed (I)=* != when travelling from the top of the atmosphere to layer i; =* != DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* != NID - INTEGER, number of layers crossed by the direct beam when (I)=* != travelling from the top of the atmosphere to layer i; =* != NID(i), i = 0..NZ-1 =* != DELTA - LOGICAL, switch to use delta-scaling (I)=* != .TRUE. -> apply delta-scaling =* != .FALSE.-> do not apply delta-scaling =* != FDR - REAL, contribution of the direct component to the total (O)=* != actinic flux at each altitude level =* != FUP - REAL, contribution of the diffuse upwelling component to (O)=* != the total actinic flux at each altitude level =* != FDN - REAL, contribution of the diffuse downwelling component to (O)=* != the total actinic flux at each altitude level =* != EDR - REAL, contribution of the direct component to the total (O)=* != spectral irradiance at each altitude level =* != EUP - REAL, contribution of the diffuse upwelling component to (O)=* != the total spectral irradiance at each altitude level =* != EDN - REAL, contribution of the diffuse downwelling component to (O)=* != the total spectral irradiance at each altitude level =* !-----------------------------------------------------------------------------* use params_mod, only : pi, precis, largest !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- INTEGER, intent(in) :: nlyr, nlev REAL, intent(in) :: zen, rsfc REAL, intent(in) :: tauu(nlyr), omu(nlyr), gu(nlyr) REAL, intent(in) :: dsdh(0:nlyr,nlyr) INTEGER, intent(in) :: nid(0:nlyr) LOGICAL, intent(in) :: delta REAL, intent(out) :: fup(nlev), fdn(nlev), fdr(nlev) REAL, intent(out) :: eup(nlev), edn(nlev), edr(nlev) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- REAL, PARAMETER :: eps = 1.E-3 !----------------------------------------------------------------------------- ! initial conditions: pi*solar flux = 1; diffuse incidence = 0 !----------------------------------------------------------------------------- REAL, PARAMETER :: pifs = 1. REAL, PARAMETER :: fdn0 = 0. REAL :: mu, sum REAL :: tausla(0:nlyr) REAL :: tauc(0:nlyr) REAL :: mu2(0:nlyr) !----------------------------------------------------------------------------- ! internal coefficients and matrix !----------------------------------------------------------------------------- INTEGER :: row, nlyrm1 REAL :: lam(nlyr), taun(nlyr), bgam(nlyr) REAL :: e1(nlyr), e2(nlyr), e3(nlyr), e4(nlyr) REAL :: cup(nlyr), cdn(nlyr), cuptn(nlyr), cdntn(nlyr) REAL :: mu1(nlyr) REAL :: a(2*nlyr), b(2*nlyr), d(2*nlyr), e(2*nlyr), y(2*nlyr) REAL :: f, g, om, tmpg REAL :: gam1, gam2, gam3, gam4 REAL :: gi(nlyr), omi(nlyr) !----------------------------------------------------------------------------- ! 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. !----------------------------------------------------------------------------- INTEGER :: mrows, mrowsm1, mrowsm2 REAL :: expon, expon0, expon1, divisr, tmp, up, dn REAL :: ssfc REAL :: wrk, wrk1 INTEGER :: i, im1, j, k !----------------------------------------------------------------------------- ! 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 ! NLYR = number of layers in the atmosphere ! NLEVEL = nlyr + 1 = number of levels !----------------------------------------------------------------------------- mu = COS( zen*pi/180. ) !----------------------------------------------------------------------------- !************* 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:nlyr) = 0. tausla(0:nlyr) = 0. mu2(0:nlyr) = 1./SQRT(largest) IF( .NOT. delta ) THEN gi(1:nlyr) = gu(1:nlyr) omi(1:nlyr) = omu(1:nlyr) taun(1:nlyr) = tauu(1:nlyr) ELSE !----------------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------------- DO k = 1, nlyr f = gu(k)*gu(k) wrk = 1. - f wrk1 = 1. - omu(k)*f gi(k) = (gu(k) - f)/wrk omi(k) = wrk*omu(k)/wrk1 taun(k) = wrk1*tauu(k) ENDDO END IF !----------------------------------------------------------------------------- ! calculate slant optical depth at the top of the atmosphere when zen>90. ! in this case, higher altitude of the top layer is recommended which can ! be easily changed in gridz.f. !----------------------------------------------------------------------------- IF( zen > 90.0 ) THEN IF(nid(0) < 0) THEN tausla(0) = largest ELSE sum = 0.0 DO j = 1, nid(0) sum = sum + 2.*taun(j)*dsdh(0,j) END DO tausla(0) = sum END IF END IF layer_loop : & DO i = 1, nlyr im1 = i - 1 g = gi(i) om = omi(i) tauc(i) = tauc(im1) + taun(i) !----------------------------------------------------------------------------- ! stay away from 1 by precision. For g, also stay away from -1 !----------------------------------------------------------------------------- tmpg = MIN( abs(g),1. - precis ) g = SIGN( tmpg,g ) om = MIN( om,1. - precis ) !----------------------------------------------------------------------------- ! calculate slant optical depth !----------------------------------------------------------------------------- IF(nid(i) < 0) THEN tausla(i) = largest ELSE sum = 0.0 DO j = 1, MIN(nid(i),i) sum = sum + taun(j)*dsdh(i,j) ENDDO DO j = MIN(nid(i),i)+1,nid(i) sum = sum + 2.*taun(j)*dsdh(i,j) ENDDO tausla(i) = sum IF(tausla(i) == tausla(im1)) THEN mu2(i) = SQRT(largest) ELSE mu2(i) = (tauc(i) - tauc(im1))/(tausla(i) - tausla(im1)) mu2(i) = SIGN( MAX( ABS(mu2(i)),1./SQRT(largest) ), mu2(i) ) END IF END IF !----------------------------------------------------------------------------- !** the following gamma equations are from pg 16,289, Table 1 !** save mu1 for each approx. for use in converting irradiance to actinic flux ! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452): !----------------------------------------------------------------------------- gam1 = .25*(7. - om*(4. + 3.*g)) gam2 = -.25*(1. - om*(4. - 3.*g)) gam3 = .25*(2. - 3.*g*mu) gam4 = 1. - gam3 mu1(i) = 0.5 !----------------------------------------------------------------------------- ! lambda = pg 16,290 equation 21 ! big gamma = pg 16,290 equation 22 ! checked limit for gam2/gam1 <<1: bgam -> (1/2)*gma2/gam1 ! so if if gam2 = 0., then bgam = 0. !----------------------------------------------------------------------------- lam(i) = sqrt(gam1*gam1 - gam2*gam2) IF( gam2 /= 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( -tausla(im1) ) expon1 = EXP( -tausla(i) ) divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i)) tmp = MAX( eps,abs(divisr) ) divisr = SIGN( tmp,divisr ) up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr dn = om*pifs*((gam1 + 1./mu2(i))*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 end do layer_loop !----------------------------------------------------------------------------- !**************** set up matrix ****** ! ssfc = pg 16,292 equation 37 where pi Fs is one (unity). !----------------------------------------------------------------------------- ssfc = rsfc*mu*EXP( -tausla(nlyr) )*pifs !----------------------------------------------------------------------------- ! MROWS = the number of rows in the matrix !----------------------------------------------------------------------------- mrows = 2*nlyr mrowsm1 = mrows - 1 mrowsm2 = mrows - 2 nlyrm1 = nlyr - 1 !----------------------------------------------------------------------------- ! the following are from pg 16,292 equations 39 - 43. ! set up first row of matrix: !----------------------------------------------------------------------------- a(1) = 0. b(1) = e1(1) d(1) = -e2(1) e(1) = fdn0 - cdn(1) !----------------------------------------------------------------------------- ! set up odd rows 3 thru (MROWS - 1): !----------------------------------------------------------------------------- a(3:mrowsm1:2) = e2(1:nlyrm1)*e3(1:nlyrm1) - e4(1:nlyrm1)*e1(1:nlyrm1) b(3:mrowsm1:2) = e1(1:nlyrm1)*e1(2:nlyr) - e3(1:nlyrm1)*e3(2:nlyr) d(3:mrowsm1:2) = e3(1:nlyrm1)*e4(2:nlyr) - e1(1:nlyrm1)*e2(2:nlyr) e(3:mrowsm1:2) = e3(1:nlyrm1)*(cup(2:nlyr) - cuptn(1:nlyrm1)) & + e1(1:nlyrm1)*(cdntn(1:nlyrm1) - cdn(2:nlyr)) !----------------------------------------------------------------------------- ! set up even rows 2 thru (MROWS - 2): !----------------------------------------------------------------------------- a(2:mrowsm2:2) = e2(2:nlyr)*e1(1:nlyrm1) - e3(1:nlyrm1)*e4(2:nlyr) b(2:mrowsm2:2) = e2(1:nlyrm1)*e2(2:nlyr) - e4(1:nlyrm1)*e4(2:nlyr) d(2:mrowsm2:2) = e1(2:nlyr)*e4(2:nlyr) - e2(2:nlyr)*e3(2:nlyr) e(2:mrowsm2:2) = (cup(2:nlyr) - cuptn(1:nlyrm1))*e2(2:nlyr) & - (cdn(2:nlyr) - cdntn(1:nlyrm1))*e4(2:nlyr) !----------------------------------------------------------------------------- ! set up last row of matrix at MROWS: !----------------------------------------------------------------------------- a(mrows) = e1(nlyr) - rsfc*e3(nlyr) b(mrows) = e2(nlyr) - rsfc*e4(nlyr) d(mrows) = 0. e(mrows) = ssfc - cuptn(nlyr) + rsfc*cdntn(nlyr) !----------------------------------------------------------------------------- ! solve tri-diagonal matrix: !----------------------------------------------------------------------------- CALL tridag( a, b, d, e, y, mrows ) !----------------------------------------------------------------------------- !*** unfold solution of matrix, compute output fluxes: !----------------------------------------------------------------------------- ! the following equations are from pg 16,291 equations 31 & 32 !----------------------------------------------------------------------------- fdr(1) = EXP( -tausla(0) ) edr(1) = mu * fdr(1) edn(1) = fdn0 eup(1) = y(1)*e3(1) - y(2)*e4(1) + cup(1) fdn(1) = edn(1)/mu1(1) fup(1) = eup(1)/mu1(1) fdr(2:nlev) = EXP( -tausla(1:nlyr) ) edr(2:nlev) = mu *fdr(2:nlev) edn(2:nlev) = y(1:mrowsm1:2)*e3(1:nlyr) + y(2:mrows:2)*e4(1:nlyr) + cdntn(1:nlyr) eup(2:nlev) = y(1:mrowsm1:2)*e1(1:nlyr) + y(2:mrows:2)*e2(1:nlyr) + cuptn(1:nlyr) fdn(2:nlev) = edn(2:nlev)/mu1(1:nlyr) fup(2:nlev) = eup(2:nlev)/mu1(1:nlyr) END SUBROUTINE ps2str SUBROUTINE tridag( a, b, c, r, u, n) !----------------------------------------------------------------------------- ! solve tridiagonal system. From Numerical Recipies, p. 40 !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- INTEGER, intent(in) :: n REAL, intent(in) :: a(n), b(n), c(n), r(n) REAL, intent(out) :: u(n) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- INTEGER :: j, jm1, jp1 REAL :: bet REAL :: gam(n) CHARACTER(len=64) :: err_msg IF (b(1) == 0.) then call wrf_error_fatal( 'tridag: zero pivot @ n == 1' ) ENDIF bet = b(1) u(1) = r(1)/bet DO j = 2, n jm1 = j - 1 gam(j) = c(jm1)/bet bet = b(j) - a(j)*gam(j) IF (bet == 0.) then write(err_msg,'(''tridag: zero pivot @ n = '',i4)') j call wrf_error_fatal( trim(err_msg) ) ENDIF u(j) = (r(j) - a(j)*u(jm1))/bet END DO DO j = n - 1, 1, -1 jp1 = j + 1 u(j) = u(j) - gam(jp1)*u(jp1) END DO END SUBROUTINE tridag END MODULE rad_trans