Radiation Scheme in CCPP
Functions/Subroutines | Variables
module_radsw_main Module Reference

This module includes ncep's modifications of the rrtm-sw radiation code from aer inc.

Functions/Subroutines

subroutine, public rswinit
 This subroutine initializes non-varying module variables, conversion factors, and look-up tables. More...
 
subroutine cldprop
 This subroutine computes the cloud optical properties for each cloudy layer and g-point interval. More...
 
subroutine mcica_subcol
 This subroutine computes the sub-colum cloud profile flag array. More...
 
subroutine setcoef
 This subroutine computes various coefficients needed in radiative transfer calculation. More...
 
subroutine spcvrtc
 This subroutine computes the shortwave radiative fluxes using two-stream method. More...
 
subroutine spcvrtm
 This subroutine computes the shortwave radiative fluxes using two-stream method of h. barder and mcica,the monte-carlo independent column approximation, for the representation of sub-grid cloud variability (i.e. cloud overlap). More...
 
subroutine swflux
 This subroutine computes the upward and downward radiation fluxes. More...
 
subroutine taumol
 
subroutine, public swrad
 This subroutine is the main sw radiation routine. More...
 

Variables

character(40), parameter vtagsw ='NCEP SW v5.1 Nov 2012 -RRTMG-SW v3.8 '
 
real(kind=kind_phys), parameter eps = 1.0e-6
 
real(kind=kind_phys), parameter oneminus = 1.0 - eps
 
real(kind=kind_phys), parameter bpade = 1.0/0.278
 
real(kind=kind_phys), parameter stpfac = 296.0/1013.0
 
real(kind=kind_phys), parameter ftiny = 1.0e-12
 
real(kind=kind_phys), parameter s0 = 1368.22
 
real(kind=kind_phys), parameter f_zero = 0.0
 
real(kind=kind_phys), parameter f_one = 1.0
 
real(kind=kind_phys), parameter amdw = con_amd/con_amw
 
real(kind=kind_phys), parameter amdo3 = con_amd/con_amo3
 
integer, dimension(nblow:nbhgh) nspa
 
integer, dimension(nblow:nbhgh) nspb
 
integer, dimension(nblow:nbhgh) idxebc
 
integer, dimension(nblow:nbhgh) idxsfc
 
integer, parameter nuvb = 27
 
logical lhswb = .false.
 
logical lhsw0 = .false.
 
logical lflxprf = .false.
 
logical lfdncmp = .false.
 
real(kind=kind_phys), dimension(0:ntbmx) exp_tbl
 
real(kind=kind_phys) heatfac
 
integer, parameter ipsdsw0 = 1
 

Function/Subroutine Documentation

subroutine module_radsw_main::cldprop ( )
private

Definition at line 1429 of file radsw_main.f.

References module_radsw_cldprtb::a0r, module_radsw_cldprtb::a0s, module_radsw_cldprtb::a1s, module_radsw_cldprtb::abari, module_radsw_cldprtb::asyice2, module_radsw_cldprtb::asyice3, module_radsw_cldprtb::asyliq1, module_radsw_cldprtb::b0r, module_radsw_cldprtb::b0s, module_radsw_cldprtb::b1s, module_radsw_cldprtb::bbari, module_radsw_cldprtb::c0r, module_radsw_cldprtb::c0s, module_radsw_cldprtb::cbari, module_radsw_cldprtb::dbari, module_radsw_cldprtb::ebari, module_radsw_cldprtb::extice2, module_radsw_cldprtb::extice3, module_radsw_cldprtb::extliq1, f_one, f_zero, module_radsw_cldprtb::fbari, ftiny, idxebc, physparam::isubcsw, physparam::iswcice, physparam::iswcliq, mcica_subcol(), module_radsw_parameters::nbdsw, module_radsw_parameters::nbhgh, module_radsw_parameters::nblow, module_radsw_parameters::ngptsw, module_radsw_cldprtb::ssaice2, module_radsw_cldprtb::ssaice3, and module_radsw_cldprtb::ssaliq1.

Referenced by swrad().

1429 ! --- inputs:
1430  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, &
1431  & cf1, nlay, ipseed, &
1432 ! --- output:
1433  & taucw, ssacw, asycw, cldfrc, cldfmc &
1434  & )
1435 
1436 ! =================== program usage description =================== !
1437 ! !
1438 ! Purpose: Compute the cloud optical properties for each cloudy layer !
1439 ! and g-point interval. !
1440 ! !
1441 ! subprograms called: none !
1442 ! !
1443 ! ==================== defination of variables ==================== !
1444 ! !
1445 ! inputs: size !
1446 ! cfrac - real, layer cloud fraction nlay !
1447 ! ..... for iswcliq > 0 (prognostic cloud sckeme) - - - !
1448 ! cliqp - real, layer in-cloud liq water path (g/m**2) nlay !
1449 ! reliq - real, mean eff radius for liq cloud (micron) nlay !
1450 ! cicep - real, layer in-cloud ice water path (g/m**2) nlay !
1451 ! reice - real, mean eff radius for ice cloud (micron) nlay !
1452 ! cdat1 - real, layer rain drop water path (g/m**2) nlay !
1453 ! cdat2 - real, effective radius for rain drop (micron) nlay !
1454 ! cdat3 - real, layer snow flake water path(g/m**2) nlay !
1455 ! cdat4 - real, mean eff radius for snow flake(micron) nlay !
1456 ! ..... for iswcliq = 0 (diagnostic cloud sckeme) - - - !
1457 ! cdat1 - real, layer cloud optical depth nlay !
1458 ! cdat2 - real, layer cloud single scattering albedo nlay !
1459 ! cdat3 - real, layer cloud asymmetry factor nlay !
1460 ! cdat4 - real, optional use nlay !
1461 ! cliqp - real, not used nlay !
1462 ! cicep - real, not used nlay !
1463 ! reliq - real, not used nlay !
1464 ! reice - real, not used nlay !
1465 ! !
1466 ! cf1 - real, effective total cloud cover at surface 1 !
1467 ! nlay - integer, vertical layer number 1 !
1468 ! ipseed- permutation seed for generating random numbers (isubcsw>0) !
1469 ! !
1470 ! outputs: !
1471 ! taucw - real, cloud optical depth, w/o delta scaled nlay*nbdsw !
1472 ! ssacw - real, weighted cloud single scattering albedo nlay*nbdsw !
1473 ! (ssa = ssacw / taucw) !
1474 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
1475 ! (asy = asycw / ssacw) !
1476 ! cldfrc - real, cloud fraction of grid mean value nlay !
1477 ! cldfmc - real, cloud fraction for each sub-column nlay*ngptsw!
1478 ! !
1479 ! !
1480 ! explanation of the method for each value of iswcliq, and iswcice. !
1481 ! set up in module "physparam" !
1482 ! !
1483 ! iswcliq=0 : input cloud optical property (tau, ssa, asy). !
1484 ! (used for diagnostic cloud method) !
1485 ! iswcliq>0 : input cloud liq/ice path and effective radius, also !
1486 ! require the user of 'iswcice' to specify the method !
1487 ! used to compute aborption due to water/ice parts. !
1488 ! ................................................................... !
1489 ! !
1490 ! iswcliq=1 : liquid water cloud optical properties are computed !
1491 ! as in hu and stamnes (1993), j. clim., 6, 728-742. !
1492 ! !
1493 ! iswcice used only when iswcliq > 0 !
1494 ! the cloud ice path (g/m2) and ice effective radius !
1495 ! (microns) are inputs. !
1496 ! iswcice=1 : ice cloud optical properties are computed as in !
1497 ! ebert and curry (1992), jgr, 97, 3831-3836. !
1498 ! iswcice=2 : ice cloud optical properties are computed as in !
1499 ! streamer v3.0 (2001), key, streamer user's guide, !
1500 ! cooperative institude for meteorological studies,95pp!
1501 ! iswcice=3 : ice cloud optical properties are computed as in !
1502 ! fu (1996), j. clim., 9. !
1503 ! !
1504 ! other cloud control module variables: !
1505 ! isubcsw =0: standard cloud scheme, no sub-col cloud approximation !
1506 ! >0: mcica sub-col cloud scheme using ipseed as permutation!
1507 ! seed for generating rundom numbers !
1508 ! !
1509 ! ====================== end of description block ================= !
1510 !
1511  use module_radsw_cldprtb
1512 
1513 ! --- inputs:
1514  integer, intent(in) :: nlay, ipseed
1515  real (kind=kind_phys), intent(in) :: cf1
1516 
1517  real (kind=kind_phys), dimension(nlay), intent(in) :: cliqp, &
1518  & reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, cfrac
1519 
1520 ! --- outputs:
1521  real (kind=kind_phys), dimension(nlay,ngptsw), intent(out) :: &
1522  & cldfmc
1523  real (kind=kind_phys), dimension(nlay,nbdsw), intent(out) :: &
1524  & taucw, ssacw, asycw
1525  real (kind=kind_phys), dimension(nlay), intent(out) :: cldfrc
1526 
1527 ! --- locals:
1528  real (kind=kind_phys), dimension(nblow:nbhgh) :: tauliq, tauice, &
1529  & ssaliq, ssaice, ssaran, ssasnw, asyliq, asyice, &
1530  & asyran, asysnw
1531  real (kind=kind_phys), dimension(nlay) :: cldf
1532 
1533  real (kind=kind_phys) :: dgeice, factor, fint, tauran, tausnw, &
1534  & cldliq, refliq, cldice, refice, cldran, cldsnw, refsnw, &
1535  & extcoliq, ssacoliq, asycoliq, extcoice, ssacoice, asycoice,&
1536  & dgesnw
1537 
1538  logical :: lcloudy(nlay,ngptsw)
1539  integer :: ia, ib, ig, jb, k, index
1540 
1541 !
1542 !===> ... begin here
1543 !
1544  do ib = 1, nbdsw
1545  do k = 1, nlay
1546  taucw(k,ib) = f_zero
1547  ssacw(k,ib) = f_one
1548  asycw(k,ib) = f_zero
1549  enddo
1550  enddo
1551 
1552 ! --- ... compute cloud radiative properties for a cloudy column
1553 
1554  lab_if_iswcliq : if (iswcliq > 0) then
1555 
1556  lab_do_k : do k = 1, nlay
1557  lab_if_cld : if (cfrac(k) > ftiny) then
1558 
1559 ! --- ... optical properties for rain and snow
1560  cldran = cdat1(k)
1561 ! refran = cdat2(k)
1562  cldsnw = cdat3(k)
1563  refsnw = cdat4(k)
1564  dgesnw = 1.0315 * refsnw ! for fu's snow formula
1565 
1566  tauran = cldran * a0r
1567 
1568 
1569 ! --- if use fu's formula it needs to be normalized by snow/ice density
1570 ! !not use snow density = 0.1 g/cm**3 = 0.1 g/(mu * m**2)
1571 ! use ice density = 0.9167 g/cm**3 = 0.9167 g/(mu * m**2)
1572 ! 1/0.9167 = 1.09087
1573 ! factor 1.5396=8/(3*sqrt(3)) converts reff to generalized ice particle size
1574 ! use newer factor value 1.0315
1575  if (cldsnw>f_zero .and. refsnw>10.0_kind_phys) then
1576 ! tausnw = cldsnw * (a0s + a1s/refsnw)
1577  tausnw = cldsnw*1.09087*(a0s + a1s/dgesnw) ! fu's formula
1578  else
1579  tausnw = f_zero
1580  endif
1581 
1582  do ib = nblow, nbhgh
1583  ssaran(ib) = tauran * (f_one - b0r(ib))
1584  ssasnw(ib) = tausnw * (f_one - (b0s(ib)+b1s(ib)*dgesnw))
1585  asyran(ib) = ssaran(ib) * c0r(ib)
1586  asysnw(ib) = ssasnw(ib) * c0s(ib)
1587  enddo
1588 
1589  cldliq = cliqp(k)
1590  cldice = cicep(k)
1591  refliq = reliq(k)
1592  refice = reice(k)
1593 
1594 ! --- ... calculation of absorption coefficients due to water clouds.
1595 
1596  if ( cldliq <= f_zero ) then
1597  do ib = nblow, nbhgh
1598  tauliq(ib) = f_zero
1599  ssaliq(ib) = f_zero
1600  asyliq(ib) = f_zero
1601  enddo
1602  else
1603  if ( iswcliq == 1 ) then
1604  factor = refliq - 1.5
1605  index = max( 1, min( 57, int( factor ) ))
1606  fint = factor - float(index)
1607 
1608  do ib = nblow, nbhgh
1609  extcoliq = max(f_zero, extliq1(index,ib) &
1610  & + fint*(extliq1(index+1,ib)-extliq1(index,ib)) )
1611  ssacoliq = max(f_zero, min(f_one, ssaliq1(index,ib) &
1612  & + fint*(ssaliq1(index+1,ib)-ssaliq1(index,ib)) ))
1613 
1614  asycoliq = max(f_zero, min(f_one, asyliq1(index,ib) &
1615  & + fint*(asyliq1(index+1,ib)-asyliq1(index,ib)) ))
1616 ! forcoliq = asycoliq * asycoliq
1617 
1618  tauliq(ib) = cldliq * extcoliq
1619  ssaliq(ib) = tauliq(ib) * ssacoliq
1620  asyliq(ib) = ssaliq(ib) * asycoliq
1621  enddo
1622  endif ! end if_iswcliq_block
1623  endif ! end if_cldliq_block
1624 
1625 ! --- ... calculation of absorption coefficients due to ice clouds.
1626 
1627  if ( cldice <= f_zero ) then
1628  do ib = nblow, nbhgh
1629  tauice(ib) = f_zero
1630  ssaice(ib) = f_zero
1631  asyice(ib) = f_zero
1632  enddo
1633  else
1634 
1635 ! --- ... ebert and curry approach for all particle sizes though somewhat
1636 ! unjustified for large ice particles
1637 
1638  if ( iswcice == 1 ) then
1639  refice = min(130.0_kind_phys,max(13.0_kind_phys,refice))
1640 
1641  do ib = nblow, nbhgh
1642  ia = idxebc(ib) ! eb_&_c band index for ice cloud coeff
1643 
1644  extcoice = max(f_zero, abari(ia)+bbari(ia)/refice )
1645  ssacoice = max(f_zero, min(f_one, &
1646  & f_one-cbari(ia)-dbari(ia)*refice ))
1647  asycoice = max(f_zero, min(f_one, &
1648  & ebari(ia)+fbari(ia)*refice ))
1649 ! forcoice = asycoice * asycoice
1650 
1651  tauice(ib) = cldice * extcoice
1652  ssaice(ib) = tauice(ib) * ssacoice
1653  asyice(ib) = ssaice(ib) * asycoice
1654  enddo
1655 
1656 ! --- ... streamer approach for ice effective radius between 5.0 and 131.0 microns
1657 
1658  elseif ( iswcice == 2 ) then
1659  refice = min(131.0_kind_phys,max(5.0_kind_phys,refice))
1660 
1661  factor = (refice - 2.0) / 3.0
1662  index = max( 1, min( 42, int( factor ) ))
1663  fint = factor - float(index)
1664 
1665  do ib = nblow, nbhgh
1666  extcoice = max(f_zero, extice2(index,ib) &
1667  & + fint*(extice2(index+1,ib)-extice2(index,ib)) )
1668  ssacoice = max(f_zero, min(f_one, ssaice2(index,ib) &
1669  & + fint*(ssaice2(index+1,ib)-ssaice2(index,ib)) ))
1670  asycoice = max(f_zero, min(f_one, asyice2(index,ib) &
1671  & + fint*(asyice2(index+1,ib)-asyice2(index,ib)) ))
1672 ! forcoice = asycoice * asycoice
1673 
1674  tauice(ib) = cldice * extcoice
1675  ssaice(ib) = tauice(ib) * ssacoice
1676  asyice(ib) = ssaice(ib) * asycoice
1677  enddo
1678 
1679 ! --- ... fu's approach for ice effective radius between 4.8 and 135 microns
1680 ! (generalized effective size from 5 to 140 microns)
1681 
1682  elseif ( iswcice == 3 ) then
1683  dgeice = max( 5.0, min( 140.0, 1.0315*refice ))
1684 
1685  factor = (dgeice - 2.0) / 3.0
1686  index = max( 1, min( 45, int( factor ) ))
1687  fint = factor - float(index)
1688 
1689  do ib = nblow, nbhgh
1690  extcoice = max(f_zero, extice3(index,ib) &
1691  & + fint*(extice3(index+1,ib)-extice3(index,ib)) )
1692  ssacoice = max(f_zero, min(f_one, ssaice3(index,ib) &
1693  & + fint*(ssaice3(index+1,ib)-ssaice3(index,ib)) ))
1694  asycoice = max(f_zero, min(f_one, asyice3(index,ib) &
1695  & + fint*(asyice3(index+1,ib)-asyice3(index,ib)) ))
1696 ! fdelta = max(f_zero, min(f_one, fdlice3(index,ib) &
1697 ! & + fint*(fdlice3(index+1,ib)-fdlice3(index,ib)) ))
1698 ! forcoice = min( asycoice, fdelta+0.5/ssacoice ) ! see fu 1996 p. 2067
1699 
1700  tauice(ib) = cldice * extcoice
1701  ssaice(ib) = tauice(ib) * ssacoice
1702  asyice(ib) = ssaice(ib) * asycoice
1703  enddo
1704 
1705  endif ! end if_iswcice_block
1706  endif ! end if_cldice_block
1707 
1708  do ib = 1, nbdsw
1709  jb = nblow + ib - 1
1710  taucw(k,ib) = tauliq(jb)+tauice(jb)+tauran+tausnw
1711  ssacw(k,ib) = ssaliq(jb)+ssaice(jb)+ssaran(jb)+ssasnw(jb)
1712  asycw(k,ib) = asyliq(jb)+asyice(jb)+asyran(jb)+asysnw(jb)
1713  enddo
1714 
1715  endif lab_if_cld
1716  enddo lab_do_k
1717 
1718  else lab_if_iswcliq
1719 
1720  do k = 1, nlay
1721  if (cfrac(k) > ftiny) then
1722  do ib = 1, nbdsw
1723  taucw(k,ib) = cdat1(k)
1724  ssacw(k,ib) = cdat1(k) * cdat2(k)
1725  asycw(k,ib) = ssacw(k,ib) * cdat3(k)
1726  enddo
1727  endif
1728  enddo
1729 
1730  endif lab_if_iswcliq
1731 
1732 ! --- distribute cloud properties to each g-point
1733 
1734  if ( isubcsw > 0 ) then ! mcica sub-col clouds approx
1735 
1736  cldf(:) = cfrac(:)
1737  where (cldf(:) < ftiny)
1738  cldf(:) = f_zero
1739  end where
1740 
1741 ! --- ... call sub-column cloud generator
1742 
1743  call mcica_subcol &
1744 ! --- inputs:
1745  & ( cldf, nlay, ipseed, &
1746 ! --- outputs:
1747  & lcloudy &
1748  & )
1749 
1750  do ig = 1, ngptsw
1751  do k = 1, nlay
1752  if ( lcloudy(k,ig) ) then
1753  cldfmc(k,ig) = f_one
1754  else
1755  cldfmc(k,ig) = f_zero
1756  endif
1757  enddo
1758  enddo
1759 
1760  else ! non-mcica, normalize cloud
1761 
1762  do k = 1, nlay
1763  cldfrc(k) = cfrac(k) / cf1
1764  enddo
1765  endif ! end if_isubcsw_block
1766 
1767  return
1768 !...................................
module_radsw_cldprtb::a1s
real(kind=kind_phys), public a1s
Definition: radsw_datatb.f:1803
module_radsw_cldprtb::ebari
real(kind=kind_phys), dimension(5), public ebari
Definition: radsw_datatb.f:186
module_radsw_cldprtb::c0r
real(kind=kind_phys), dimension(nblow:nbhgh), public c0r
Definition: radsw_datatb.f:1806
module_radsw_cldprtb::asyice2
real(kind=kind_phys), dimension(43, nblow:nbhgh), public asyice2
Definition: radsw_datatb.f:182
module_radsw_cldprtb::extice2
real(kind=kind_phys), dimension(43, nblow:nbhgh), public extice2
Definition: radsw_datatb.f:182
module_radsw_cldprtb::a0s
real(kind=kind_phys), public a0s
Definition: radsw_datatb.f:1803
module_radsw_cldprtb::extliq1
real(kind=kind_phys), dimension(58, nblow:nbhgh), public extliq1
Definition: radsw_datatb.f:180
module_radsw_cldprtb
This module contains coefficients of cloud optical properties for each of the spectral bands...
Definition: radsw_datatb.f:126
module_radsw_cldprtb::c0s
real(kind=kind_phys), dimension(nblow:nbhgh), public c0s
Definition: radsw_datatb.f:1806
module_radsw_cldprtb::b0s
real(kind=kind_phys), dimension(nblow:nbhgh), public b0s
Definition: radsw_datatb.f:1806
module_radsw_cldprtb::ssaice2
real(kind=kind_phys), dimension(43, nblow:nbhgh), public ssaice2
Definition: radsw_datatb.f:182
module_radsw_cldprtb::extice3
real(kind=kind_phys), dimension(46, nblow:nbhgh), public extice3
Definition: radsw_datatb.f:184
module_radsw_cldprtb::fbari
real(kind=kind_phys), dimension(5), public fbari
Definition: radsw_datatb.f:186
module_radsw_cldprtb::dbari
real(kind=kind_phys), dimension(5), public dbari
Definition: radsw_datatb.f:186
module_radsw_cldprtb::asyliq1
real(kind=kind_phys), dimension(58, nblow:nbhgh), public asyliq1
Definition: radsw_datatb.f:180
module_radsw_cldprtb::asyice3
real(kind=kind_phys), dimension(46, nblow:nbhgh), public asyice3
Definition: radsw_datatb.f:184
module_radsw_cldprtb::ssaliq1
real(kind=kind_phys), dimension(58, nblow:nbhgh), public ssaliq1
Definition: radsw_datatb.f:180
module_radsw_cldprtb::ssaice3
real(kind=kind_phys), dimension(46, nblow:nbhgh), public ssaice3
Definition: radsw_datatb.f:184
module_radsw_cldprtb::a0r
real(kind=kind_phys), public a0r
Definition: radsw_datatb.f:1803
module_radsw_cldprtb::cbari
real(kind=kind_phys), dimension(5), public cbari
Definition: radsw_datatb.f:186
module_radsw_cldprtb::b1s
real(kind=kind_phys), dimension(nblow:nbhgh), public b1s
Definition: radsw_datatb.f:1806
module_radsw_cldprtb::b0r
real(kind=kind_phys), dimension(nblow:nbhgh), public b0r
Definition: radsw_datatb.f:1806
module_radsw_cldprtb::abari
real(kind=kind_phys), dimension(5), public abari
Definition: radsw_datatb.f:186
module_radsw_cldprtb::bbari
real(kind=kind_phys), dimension(5), public bbari
Definition: radsw_datatb.f:186

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine module_radsw_main::mcica_subcol ( )
Parameters
[in]cldfreal, (nlay), layer cloud fraction
[in]nlayinteger, 1, number of model vertical layers
[in]ipseedinteger, 1, permute seed for random num generator
[out]lcloudylogical, (nlay,ngptsw), sub-colum cloud profile flag array

External Module Variables


physparam::iovrsw - control flag for cloud overlapping method
=0:random; =1:maximum/random; =2:maximum

Definition at line 1783 of file radsw_main.f.

References f_one, physparam::iovrsw, and module_radsw_parameters::ngptsw.

Referenced by cldprop().

1783 ! --- inputs:
1784  & ( cldf, nlay, ipseed, &
1785 ! --- outputs:
1786  & lcloudy &
1787  & )
1788 
1789 ! ==================== defination of variables ==================== !
1790 ! !
1791 ! input variables: size !
1792 ! cldf - real, layer cloud fraction nlay !
1793 ! nlay - integer, number of model vertical layers 1 !
1794 ! ipseed - integer, permute seed for random num generator 1 !
1795 ! ** note : if the cloud generator is called multiple times, need !
1796 ! to permute the seed between each call; if between calls !
1797 ! for lw and sw, use values differ by the number of g-pts. !
1798 ! !
1799 ! output variables: !
1800 ! lcloudy - logical, sub-colum cloud profile flag array nlay*ngptsw!
1801 ! !
1802 ! other control flags from module variables: !
1803 ! iovrsw : control flag for cloud overlapping method !
1804 ! =0:random; =1:maximum/random; =2:maximum !
1805 ! !
1806 ! !
1807 ! ===================== end of definitions ==================== !
1808 
1809  implicit none
1810 
1811 ! --- inputs:
1812  integer, intent(in) :: nlay, ipseed
1813 
1814  real (kind=kind_phys), dimension(nlay), intent(in) :: cldf
1815 
1816 ! --- outputs:
1817  logical, dimension(nlay,ngptsw), intent(out):: lcloudy
1818 
1819 ! --- locals:
1820  real (kind=kind_phys) :: cdfunc(nlay,ngptsw), tem1, &
1821  & rand2d(nlay*ngptsw), rand1d(ngptsw)
1822 
1823  type(random_stat) :: stat ! for thread safe random generator
1824 
1825  integer :: k, n, k1
1826 !
1827 !===> ... begin here
1828 !
1829 ! --- ... advance randum number generator by ipseed values
1830 
1831  call random_setseed &
1832 ! --- inputs:
1833  & ( ipseed, &
1834 ! --- outputs:
1835  & stat &
1836  & )
1837 
1838 ! --- ... sub-column set up according to overlapping assumption
1839 
1840  select case ( iovrsw )
1841 
1842  case( 0 ) ! random overlap, pick a random value at every level
1843 
1844  call random_number &
1845 ! --- inputs: ( none )
1846 ! --- outputs:
1847  & ( rand2d, stat )
1848 
1849  k1 = 0
1850  do n = 1, ngptsw
1851  do k = 1, nlay
1852  k1 = k1 + 1
1853  cdfunc(k,n) = rand2d(k1)
1854  enddo
1855  enddo
1856 
1857  case( 1 ) ! max-ran overlap
1858 
1859  call random_number &
1860 ! --- inputs: ( none )
1861 ! --- outputs:
1862  & ( rand2d, stat )
1863 
1864  k1 = 0
1865  do n = 1, ngptsw
1866  do k = 1, nlay
1867  k1 = k1 + 1
1868  cdfunc(k,n) = rand2d(k1)
1869  enddo
1870  enddo
1871 
1872 ! --- first pick a random number for bottom/top layer.
1873 ! then walk up the column: (aer's code)
1874 ! if layer below is cloudy, use the same rand num in the layer below
1875 ! if layer below is clear, use a new random number
1876 
1877 ! --- from bottom up
1878  do k = 2, nlay
1879  k1 = k - 1
1880  tem1 = f_one - cldf(k1)
1881 
1882  do n = 1, ngptsw
1883  if ( cdfunc(k1,n) > tem1 ) then
1884  cdfunc(k,n) = cdfunc(k1,n)
1885  else
1886  cdfunc(k,n) = cdfunc(k,n) * tem1
1887  endif
1888  enddo
1889  enddo
1890 
1891 ! --- then walk down the column: (if use original author's method)
1892 ! if layer above is cloudy, use the same rand num in the layer above
1893 ! if layer above is clear, use a new random number
1894 
1895 ! --- from top down
1896 ! do k = nlay-1, 1, -1
1897 ! k1 = k + 1
1898 ! tem1 = f_one - cldf(k1)
1899 
1900 ! do n = 1, ngptsw
1901 ! if ( cdfunc(k1,n) > tem1 ) then
1902 ! cdfunc(k,n) = cdfunc(k1,n)
1903 ! else
1904 ! cdfunc(k,n) = cdfunc(k,n) * tem1
1905 ! endif
1906 ! enddo
1907 ! enddo
1908 
1909  case( 2 ) ! maximum overlap, pick same random numebr at every level
1910 
1911  call random_number &
1912 ! --- inputs: ( none )
1913 ! --- outputs:
1914  & ( rand1d, stat )
1915 
1916  do n = 1, ngptsw
1917  tem1 = rand1d(n)
1918 
1919  do k = 1, nlay
1920  cdfunc(k,n) = tem1
1921  enddo
1922  enddo
1923 
1924  end select
1925 
1926 ! --- ... generate subcolumns for homogeneous clouds
1927 
1928  do k = 1, nlay
1929  tem1 = f_one - cldf(k)
1930 
1931  do n = 1, ngptsw
1932  lcloudy(k,n) = cdfunc(k,n) >= tem1
1933  enddo
1934  enddo
1935 
1936  return
1937 ! ..................................

Here is the caller graph for this function:

subroutine, public module_radsw_main::rswinit ( )
Parameters
[in]meinteger, print control for parallel process
[out]NONE

External Module Variables


physparam::iswrate - heating rate unit selections
=1: output in k/day
=2: output in k/second
physparam::iswrgas - control flag for rare gases (ch4,n2o,o2, etc.)
=0: do not include rare gases
>0: include all rare gases
physparam::iswcliq - liquid cloud optical properties contrl flag
=0: input cloud opt depth from diagnostic scheme
>0: input cwp,rew, and other cloud content parameters
physparam::isubcsw - sub-column cloud approximation control flag
=0: no sub-col cld treatment, use grid-mean cld quantities
=1: mcica sub-col, prescribed seeds to get random numbers
=2: mcica sub-col, providing array icseed for random numbers
physparam::icldflg - cloud scheme control flag
=0: diagnostic scheme gives cloud tau, omiga, and g.
=1: prognostic scheme gives cloud liq/ice path, etc.
physparam::iovrsw - clouds vertical overlapping control flag
=0: random overlapping clouds
=1: maximum/random overlapping clouds
=2: maximum overlap cloud
physparam::iswmode - control flag for 2-stream transfer scheme
=1; delta-eddington (joseph et al., 1976)
=2: pifm (zdunkowski et al., 1980)
=3: discrete ordinates (liou, 1973)

Definition at line 1272 of file radsw_main.f.

References bpade, physcons::con_cp, physcons::con_g, exp_tbl, heatfac, physparam::icldflg, physparam::iovrsw, physparam::isubcsw, physparam::iswcliq, physparam::iswmode, physparam::iswrate, physparam::iswrgas, module_radsw_parameters::ntbmx, and vtagsw.

Referenced by module_radiation_driver::radinit().

1272 
1273 ! --- inputs:
1274  & ( me )
1275 ! --- outputs: (none)
1276 
1277 ! =================== program usage description =================== !
1278 ! !
1279 ! purpose: initialize non-varying module variables, conversion factors,!
1280 ! and look-up tables. !
1281 ! !
1282 ! subprograms called: none !
1283 ! !
1284 ! ==================== defination of variables ==================== !
1285 ! !
1286 ! inputs: !
1287 ! me - print control for parallel process !
1288 ! !
1289 ! outputs: (none) !
1290 ! !
1291 ! external module variables: (in physparam) !
1292 ! iswrate - heating rate unit selections !
1293 ! =1: output in k/day !
1294 ! =2: output in k/second !
1295 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
1296 ! =0: do not include rare gases !
1297 ! >0: include all rare gases !
1298 ! iswcliq - liquid cloud optical properties contrl flag !
1299 ! =0: input cloud opt depth from diagnostic scheme !
1300 ! >0: input cwp,rew, and other cloud content parameters !
1301 ! isubcsw - sub-column cloud approximation control flag !
1302 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
1303 ! =1: mcica sub-col, prescribed seeds to get random numbers !
1304 ! =2: mcica sub-col, providing array icseed for random numbers!
1305 ! icldflg - cloud scheme control flag !
1306 ! =0: diagnostic scheme gives cloud tau, omiga, and g. !
1307 ! =1: prognostic scheme gives cloud liq/ice path, etc. !
1308 ! iovrsw - clouds vertical overlapping control flag !
1309 ! =0: random overlapping clouds !
1310 ! =1: maximum/random overlapping clouds !
1311 ! =2: maximum overlap cloud !
1312 ! iswmode - control flag for 2-stream transfer scheme !
1313 ! =1; delta-eddington (joseph et al., 1976) !
1314 ! =2: pifm (zdunkowski et al., 1980) !
1315 ! =3: discrete ordinates (liou, 1973) !
1316 ! !
1317 ! ******************************************************************* !
1318 ! !
1319 ! definitions: !
1320 ! arrays for 10000-point look-up tables: !
1321 ! tau_tbl clear-sky optical depth !
1322 ! exp_tbl exponential lookup table for transmittance !
1323 ! !
1324 ! ******************************************************************* !
1325 ! !
1326 ! ====================== end of description block ================= !
1327 
1328 ! --- inputs:
1329  integer, intent(in) :: me
1330 
1331 ! --- outputs: none
1332 
1333 ! --- locals:
1334  real (kind=kind_phys), parameter :: expeps = 1.e-20
1335 
1336  integer :: i
1337 
1338  real (kind=kind_phys) :: tfn, tau
1339 
1340 !
1341 !===> ... begin here
1342 !
1343  if ( iovrsw<0 .or. iovrsw>2 ) then
1344  print *,' *** Error in specification of cloud overlap flag', &
1345  & ' IOVRSW=',iovrsw,' in RSWINIT !!'
1346  stop
1347  endif
1348 
1349  if (me == 0) then
1350  print *,' - Using AER Shortwave Radiation, Version: ',vtagsw
1351 
1352  if (iswmode == 1) then
1353  print *,' --- Delta-eddington 2-stream transfer scheme'
1354  else if (iswmode == 2) then
1355  print *,' --- PIFM 2-stream transfer scheme'
1356  else if (iswmode == 3) then
1357  print *,' --- Discrete ordinates 2-stream transfer scheme'
1358  endif
1359 
1360  if (iswrgas <= 0) then
1361  print *,' --- Rare gases absorption is NOT included in SW'
1362  else
1363  print *,' --- Include rare gases N2O, CH4, O2, absorptions',&
1364  & ' in SW'
1365  endif
1366 
1367  if ( isubcsw == 0 ) then
1368  print *,' --- Using standard grid average clouds, no ', &
1369  & 'sub-column clouds approximation applied'
1370  elseif ( isubcsw == 1 ) then
1371  print *,' --- Using MCICA sub-colum clouds approximation ', &
1372  & 'with a prescribed sequence of permutation seeds'
1373  elseif ( isubcsw == 2 ) then
1374  print *,' --- Using MCICA sub-colum clouds approximation ', &
1375  & 'with provided input array of permutation seeds'
1376  else
1377  print *,' *** Error in specification of sub-column cloud ', &
1378  & ' control flag isubcsw =',isubcsw,' !!'
1379  stop
1380  endif
1381  endif
1382 
1383 ! --- ... check cloud flags for consistency
1384 
1385  if ((icldflg == 0 .and. iswcliq /= 0) .or. &
1386  & (icldflg == 1 .and. iswcliq == 0)) then
1387  print *,' *** Model cloud scheme inconsistent with SW', &
1388  & ' radiation cloud radiative property setup !!'
1389  stop
1390  endif
1391 
1392 ! --- ... setup constant factors for heating rate
1393 ! the 1.0e-2 is to convert pressure from mb to N/m**2
1394 
1395  if (iswrate == 1) then
1396 ! heatfac = 8.4391
1397 ! heatfac = con_g * 86400. * 1.0e-2 / con_cp ! (in k/day)
1398  heatfac = con_g * 864.0 / con_cp ! (in k/day)
1399  else
1400  heatfac = con_g * 1.0e-2 / con_cp ! (in k/second)
1401  endif
1402 
1403 ! --- ... define exponential lookup tables for transmittance. tau is
1404 ! computed as a function of the tau transition function, and
1405 ! transmittance is calculated as a function of tau. all tables
1406 ! are computed at intervals of 0.0001. the inverse of the
1407 ! constant used in the Pade approximation to the tau transition
1408 ! function is set to bpade.
1409 
1410  exp_tbl(0) = 1.0
1411  exp_tbl(ntbmx) = expeps
1412 
1413  do i = 1, ntbmx-1
1414  tfn = float(i) / float(ntbmx-i)
1415  tau = bpade * tfn
1416  exp_tbl(i) = exp( -tau )
1417  enddo
1418 
1419  return
1420 !...................................

Here is the caller graph for this function:

subroutine module_radsw_main::setcoef ( )
private

Definition at line 1945 of file radsw_main.f.

References f_one, f_zero, module_radsw_ref::preflog, stpfac, and module_radsw_ref::tref.

Referenced by swrad().

1945 ! --- inputs:
1946  & ( pavel,tavel,h2ovmr, nlay,nlp1, &
1947 ! --- outputs:
1948  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, &
1949  & selffac,selffrac,indself,forfac,forfrac,indfor &
1950  & )
1951 
1952 ! =================== program usage description =================== !
1953 ! !
1954 ! purpose: compute various coefficients needed in radiative transfer !
1955 ! calculations. !
1956 ! !
1957 ! subprograms called: none !
1958 ! !
1959 ! ==================== defination of variables ==================== !
1960 ! !
1961 ! inputs: -size- !
1962 ! pavel - real, layer pressures (mb) nlay !
1963 ! tavel - real, layer temperatures (k) nlay !
1964 ! h2ovmr - real, layer w.v. volum mixing ratio (kg/kg) nlay !
1965 ! nlay/nlp1 - integer, total number of vertical layers, levels 1 !
1966 ! !
1967 ! outputs: !
1968 ! laytrop - integer, tropopause layer index (unitless) 1 !
1969 ! jp - real, indices of lower reference pressure nlay !
1970 ! jt, jt1 - real, indices of lower reference temperatures nlay !
1971 ! at levels of jp and jp+1 !
1972 ! facij - real, factors multiply the reference ks, nlay !
1973 ! i,j=0/1 for lower/higher of the 2 appropriate !
1974 ! temperatures and altitudes. !
1975 ! selffac - real, scale factor for w. v. self-continuum nlay !
1976 ! equals (w. v. density)/(atmospheric density !
1977 ! at 296k and 1013 mb) !
1978 ! selffrac - real, factor for temperature interpolation of nlay !
1979 ! reference w. v. self-continuum data !
1980 ! indself - integer, index of lower ref temp for selffac nlay !
1981 ! forfac - real, scale factor for w. v. foreign-continuum nlay !
1982 ! forfrac - real, factor for temperature interpolation of nlay !
1983 ! reference w.v. foreign-continuum data !
1984 ! indfor - integer, index of lower ref temp for forfac nlay !
1985 ! !
1986 ! ====================== end of definitions =================== !
1987 
1988 ! --- inputs:
1989  integer, intent(in) :: nlay, nlp1
1990 
1991  real (kind=kind_phys), dimension(:), intent(in) :: pavel, tavel, &
1992  & h2ovmr
1993 
1994 ! --- outputs:
1995  integer, dimension(nlay), intent(out) :: indself, indfor, &
1996  & jp, jt, jt1
1997  integer, intent(out) :: laytrop
1998 
1999  real (kind=kind_phys), dimension(nlay), intent(out) :: fac00, &
2000  & fac01, fac10, fac11, selffac, selffrac, forfac, forfrac
2001 
2002 ! --- locals:
2003  real (kind=kind_phys) :: plog, fp, fp1, ft, ft1, tem1, tem2
2004 
2005  integer :: i, k, jp1
2006 !
2007 !===> ... begin here
2008 !
2009  laytrop= nlay
2010 
2011  do k = 1, nlay
2012 
2013  forfac(k) = pavel(k)*stpfac / (tavel(k)*(f_one + h2ovmr(k)))
2014 
2015 ! --- ... find the two reference pressures on either side of the
2016 ! layer pressure. store them in jp and jp1. store in fp the
2017 ! fraction of the difference (in ln(pressure)) between these
2018 ! two values that the layer pressure lies.
2019 
2020  plog = log(pavel(k))
2021  jp(k) = max(1, min(58, int(36.0 - 5.0*(plog+0.04)) ))
2022  jp1 = jp(k) + 1
2023  fp = 5.0 * (preflog(jp(k)) - plog)
2024 
2025 ! --- ... determine, for each reference pressure (jp and jp1), which
2026 ! reference temperature (these are different for each reference
2027 ! pressure) is nearest the layer temperature but does not exceed it.
2028 ! store these indices in jt and jt1, resp. store in ft (resp. ft1)
2029 ! the fraction of the way between jt (jt1) and the next highest
2030 ! reference temperature that the layer temperature falls.
2031 
2032  tem1 = (tavel(k) - tref(jp(k))) / 15.0
2033  tem2 = (tavel(k) - tref(jp1 )) / 15.0
2034  jt(k) = max(1, min(4, int(3.0 + tem1) ))
2035  jt1(k) = max(1, min(4, int(3.0 + tem2) ))
2036  ft = tem1 - float(jt(k) - 3)
2037  ft1 = tem2 - float(jt1(k) - 3)
2038 
2039 ! --- ... we have now isolated the layer ln pressure and temperature,
2040 ! between two reference pressures and two reference temperatures
2041 ! (for each reference pressure). we multiply the pressure
2042 ! fraction fp with the appropriate temperature fractions to get
2043 ! the factors that will be needed for the interpolation that yields
2044 ! the optical depths (performed in routines taugbn for band n).
2045 
2046  fp1 = f_one - fp
2047  fac10(k) = fp1 * ft
2048  fac00(k) = fp1 * (f_one - ft)
2049  fac11(k) = fp * ft1
2050  fac01(k) = fp * (f_one - ft1)
2051 
2052 ! --- ... if the pressure is less than ~100mb, perform a different
2053 ! set of species interpolations.
2054 
2055  if ( plog > 4.56 ) then
2056 
2057  laytrop = k
2058 
2059 ! --- ... set up factors needed to separately include the water vapor
2060 ! foreign-continuum in the calculation of absorption coefficient.
2061 
2062  tem1 = (332.0 - tavel(k)) / 36.0
2063  indfor(k) = min(2, max(1, int(tem1)))
2064  forfrac(k) = tem1 - float(indfor(k))
2065 
2066 ! --- ... set up factors needed to separately include the water vapor
2067 ! self-continuum in the calculation of absorption coefficient.
2068 
2069  tem2 = (tavel(k) - 188.0) / 7.2
2070  indself(k) = min(9, max(1, int(tem2)-7))
2071  selffrac(k) = tem2 - float(indself(k) + 7)
2072  selffac(k) = h2ovmr(k) * forfac(k)
2073 
2074  else
2075 
2076 ! --- ... set up factors needed to separately include the water vapor
2077 ! foreign-continuum in the calculation of absorption coefficient.
2078 
2079  tem1 = (tavel(k) - 188.0) / 36.0
2080  indfor(k) = 3
2081  forfrac(k) = tem1 - f_one
2082 
2083  indself(k) = 0
2084  selffrac(k) = f_zero
2085  selffac(k) = f_zero
2086 
2087  endif
2088 
2089  enddo ! end_do_k_loop
2090 
2091  return
2092 ! ..................................

Here is the caller graph for this function:

subroutine module_radsw_main::spcvrtc ( )
private

Definition at line 2100 of file radsw_main.f.

References bpade, eps, exp_tbl, f_one, f_zero, ftiny, idxsfc, physparam::iswmode, module_radsw_parameters::nbdsw, module_radsw_parameters::nblow, module_radsw_parameters::ngb, module_radsw_parameters::ngptsw, module_radsw_parameters::ntbmx, nuvb, oneminus, and swflux().

Referenced by swrad().

2100 ! --- inputs:
2101  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfrc, &
2102  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
2103  & nlay, nlp1, &
2104 ! --- outputs:
2105  & fxupc,fxdnc,fxup0,fxdn0, &
2106  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
2107  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
2108  & )
2109 
2110 ! =================== program usage description =================== !
2111 ! !
2112 ! purpose: computes the shortwave radiative fluxes using two-stream !
2113 ! method !
2114 ! !
2115 ! subprograms called: swflux !
2116 ! !
2117 ! ==================== defination of variables ==================== !
2118 ! !
2119 ! inputs: size !
2120 ! ssolar - real, incoming solar flux at top 1 !
2121 ! cosz - real, cosine solar zenith angle 1 !
2122 ! sntz - real, secant solar zenith angle 1 !
2123 ! albbm - real, surface albedo for direct beam radiation 2 !
2124 ! albdf - real, surface albedo for diffused radiation 2 !
2125 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
2126 ! cldfrc - real, layer cloud fraction nlay !
2127 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
2128 ! cf0 - real, =1-cf1 1 !
2129 ! taug - real, spectral optical depth for gases nlay*ngptsw!
2130 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
2131 ! tauae - real, aerosols optical depth nlay*nbdsw !
2132 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
2133 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
2134 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
2135 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
2136 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
2137 ! nlay,nlp1 - integer, number of layers/levels 1 !
2138 ! !
2139 ! output variables: !
2140 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
2141 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
2142 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
2143 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
2144 ! ftoauc - real, tot sky toa upwd flux 1 !
2145 ! ftoau0 - real, clr sky toa upwd flux 1 !
2146 ! ftoadc - real, toa downward (incoming) solar flux 1 !
2147 ! fsfcuc - real, tot sky sfc upwd flux 1 !
2148 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
2149 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
2150 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
2151 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
2152 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
2153 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
2154 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
2155 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
2156 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
2157 ! !
2158 ! internal variables: !
2159 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
2160 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
2161 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
2162 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
2163 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
2164 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
2165 ! !
2166 ! control parameters in module "physparam" !
2167 ! iswmode - control flag for 2-stream transfer schemes !
2168 ! = 1 delta-eddington (joseph et al., 1976) !
2169 ! = 2 pifm (zdunkowski et al., 1980) !
2170 ! = 3 discrete ordinates (liou, 1973) !
2171 ! !
2172 ! ******************************************************************* !
2173 ! original code description !
2174 ! !
2175 ! method: !
2176 ! ------- !
2177 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
2178 ! kmodts = 1 eddington (joseph et al., 1976) !
2179 ! = 2 pifm (zdunkowski et al., 1980) !
2180 ! = 3 discrete ordinates (liou, 1973) !
2181 ! !
2182 ! modifications: !
2183 ! -------------- !
2184 ! original: h. barker !
2185 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
2186 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
2187 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
2188 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
2189 ! revision: code modified so that delta scaling is not done in cloudy !
2190 ! profiles if routine cldprop is used; delta scaling can be !
2191 ! applied by swithcing code below if cldprop is not used to !
2192 ! get cloud properties. aer, jan 2005 !
2193 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
2194 ! revision: use exponential lookup table for transmittance: mjiacono, !
2195 ! aer, aug 2007 !
2196 ! !
2197 ! ******************************************************************* !
2198 ! ====================== end of description block ================= !
2199 
2200 ! --- constant parameters:
2201  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
2202  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
2203  real (kind=kind_phys), parameter :: od_lo = 0.06
2204  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
2205 
2206 ! --- inputs:
2207  integer, intent(in) :: nlay, nlp1
2208 
2209  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
2210  & taug, taur
2211  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
2212  & taucw, ssacw, asycw, tauae, ssaae, asyae
2213 
2214  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
2215  real (kind=kind_phys), dimension(nlay), intent(in) :: cldfrc
2216 
2217  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
2218 
2219  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
2220 
2221 ! --- outputs:
2222  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
2223  & fxupc, fxdnc, fxup0, fxdn0
2224 
2225  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
2226  & sfbm0, sfdf0
2227 
2228  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
2229  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
2230 
2231 ! --- locals:
2232  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
2233  & zldbt0
2234 
2235  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
2236  & ztrad, ztdbt, zldbt, zfu, zfd
2237 
2238  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
2239  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
2240  & zc0, zc1, za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, &
2241  & zrpp, zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, &
2242  & zexp3, zexp4, zden1, ze1r45, ftind, zsolar, zrefb1, &
2243  & zrefd1, ztrab1, ztrad1, ztdbt0, zr1, zr2, zr3, zr4, zr5, &
2244  & zt1, zt2, zt3, zf1, zf2
2245 
2246  integer :: ib, ibd, jb, jg, k, kp, itind
2247 !
2248 !===> ... begin here
2249 !
2250 ! --- ... initialization of output fluxes
2251 
2252  do ib = 1, nbdsw
2253  do k = 1, nlp1
2254  fxdnc(k,ib) = f_zero
2255  fxupc(k,ib) = f_zero
2256  fxdn0(k,ib) = f_zero
2257  fxup0(k,ib) = f_zero
2258  enddo
2259  enddo
2260 
2261  ftoadc = f_zero
2262  ftoauc = f_zero
2263  ftoau0 = f_zero
2264  fsfcuc = f_zero
2265  fsfcu0 = f_zero
2266  fsfcdc = f_zero
2267  fsfcd0 = f_zero
2268 
2269 !! --- ... uv-b surface downward fluxes
2270  suvbfc = f_zero
2271  suvbf0 = f_zero
2272 
2273 !! --- ... output surface flux components
2274  sfbmc(1) = f_zero
2275  sfbmc(2) = f_zero
2276  sfdfc(1) = f_zero
2277  sfdfc(2) = f_zero
2278  sfbm0(1) = f_zero
2279  sfbm0(2) = f_zero
2280  sfdf0(1) = f_zero
2281  sfdf0(2) = f_zero
2282 
2283 ! --- ... loop over all g-points in each band
2284 
2285  lab_do_jg : do jg = 1, ngptsw
2286 
2287  jb = ngb(jg)
2288  ib = jb + 1 - nblow
2289  ibd = idxsfc(jb)
2290 
2291  zsolar = ssolar * sfluxzen(jg)
2292 
2293 ! --- ... set up toa direct beam and surface values (beam and diff)
2294 
2295  ztdbt(nlp1) = f_one
2296  ztdbt0 = f_one
2297 
2298  zldbt(1) = f_zero
2299  if (ibd /= 0) then
2300  zrefb(1) = albbm(ibd)
2301  zrefd(1) = albdf(ibd)
2302  else
2303  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
2304  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
2305  endif
2306  ztrab(1) = f_zero
2307  ztrad(1) = f_zero
2308 
2309 ! --- ... compute clear-sky optical parameters, layer reflectance and transmittance
2310 
2311  do k = nlay, 1, -1
2312  kp = k + 1
2313 
2314  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
2315  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
2316  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
2317  zssaw = min( oneminus, zssa0 / ztau0 )
2318  zasyw = zasy0 / max( ftiny, zssa0 )
2319 
2320 ! --- ... saving clear-sky quantities for later total-sky usage
2321  ztaus(k) = ztau0
2322  zssas(k) = zssa0
2323  zasys(k) = zasy0
2324 
2325 ! --- ... delta scaling for clear-sky condition
2326  za1 = zasyw * zasyw
2327  za2 = zssaw * za1
2328 
2329  ztau1 = (f_one - za2) * ztau0
2330  zssa1 = (zssaw - za2) / (f_one - za2)
2331 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
2332  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
2333  zasy3 = 0.75 * zasy1
2334 
2335 ! --- ... general two-stream expressions
2336  if ( iswmode == 1 ) then
2337  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2338  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
2339  zgam3 = 0.5 - zasy3 * cosz
2340  elseif ( iswmode == 2 ) then ! pifm
2341  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2342  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2343  zgam3 = 0.5 - zasy3 * cosz
2344  elseif ( iswmode == 3 ) then ! discrete ordinates
2345  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2346  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2347  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2348  endif
2349  zgam4 = f_one - zgam3
2350 
2351 ! --- ... compute homogeneous reflectance and transmittance
2352 
2353  if ( zssaw >= zcrit ) then ! for conservative scattering
2354  za1 = zgam1 * cosz - zgam3
2355  za2 = zgam1 * ztau1
2356 
2357 ! --- ... use exponential lookup table for transmittance, or expansion
2358 ! of exponential for low optical depth
2359 
2360  zb1 = min( ztau1*sntz , 500.0 )
2361  if ( zb1 <= od_lo ) then
2362  zb2 = f_one - zb1 + 0.5*zb1*zb1
2363  else
2364  ftind = zb1 / (bpade + zb1)
2365  itind = ftind*ntbmx + 0.5
2366  zb2 = exp_tbl(itind)
2367  endif
2368 
2369 ! ... collimated beam
2370  zrefb(kp) = max(f_zero, min(f_one, &
2371  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2372  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
2373 
2374 ! ... isotropic incidence
2375  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
2376  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
2377 
2378  else ! for non-conservative scattering
2379  za1 = zgam1*zgam4 + zgam2*zgam3
2380  za2 = zgam1*zgam3 + zgam2*zgam4
2381  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2382  zrk2= 2.0 * zrk
2383 
2384  zrp = zrk * cosz
2385  zrp1 = f_one + zrp
2386  zrm1 = f_one - zrp
2387  zrpp = f_one - zrp*zrp
2388  zrkg1= zrk + zgam1
2389  zrkg3= zrk * zgam3
2390  zrkg4= zrk * zgam4
2391 
2392  zr1 = zrm1 * (za2 + zrkg3)
2393  zr2 = zrp1 * (za2 - zrkg3)
2394  zr3 = zrk2 * (zgam3 - za2*cosz)
2395  zr4 = zrpp * zrkg1
2396  zr5 = zrpp * (zrk - zgam1)
2397 
2398  zt1 = zrp1 * (za1 + zrkg4)
2399  zt2 = zrm1 * (za1 - zrkg4)
2400  zt3 = zrk2 * (zgam4 + za1*cosz)
2401 
2402 ! --- ... use exponential lookup table for transmittance, or expansion
2403 ! of exponential for low optical depth
2404 
2405  zb1 = min( zrk*ztau1, 500.0 )
2406  if ( zb1 <= od_lo ) then
2407  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2408  else
2409  ftind = zb1 / (bpade + zb1)
2410  itind = ftind*ntbmx + 0.5
2411  zexm1 = exp_tbl(itind)
2412  endif
2413  zexp1 = f_one / zexm1
2414 
2415  zb2 = min( sntz*ztau1, 500.0 )
2416  if ( zb2 <= od_lo ) then
2417  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2418  else
2419  ftind = zb2 / (bpade + zb2)
2420  itind = ftind*ntbmx + 0.5
2421  zexm2 = exp_tbl(itind)
2422  endif
2423  zexp2 = f_one / zexm2
2424  ze1r45 = zr4*zexp1 + zr5*zexm1
2425 
2426 ! ... collimated beam
2427  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
2428  zrefb(kp) = eps1
2429  ztrab(kp) = zexm2
2430  else
2431  zden1 = zssa1 / ze1r45
2432  zrefb(kp) = max(f_zero, min(f_one, &
2433  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
2434  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
2435  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
2436  endif
2437 
2438 ! ... diffuse beam
2439  zden1 = zr4 / (ze1r45 * zrkg1)
2440  zrefd(kp) = max(f_zero, min(f_one, &
2441  & zgam2*(zexp1 - zexm1)*zden1 ))
2442  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2443  endif ! end if_zssaw_block
2444 
2445 ! --- ... direct beam transmittance. use exponential lookup table
2446 ! for transmittance, or expansion of exponential for low
2447 ! optical depth
2448 
2449  zr1 = ztau1 * sntz
2450  if ( zr1 <= od_lo ) then
2451  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2452  else
2453  ftind = zr1 / (bpade + zr1)
2454  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2455  zexp3 = exp_tbl(itind)
2456  endif
2457 
2458  ztdbt(k) = zexp3 * ztdbt(kp)
2459  zldbt(kp) = zexp3
2460 
2461 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2462 ! (must use 'orig', unscaled cloud optical depth)
2463 
2464  zr1 = ztau0 * sntz
2465  if ( zr1 <= od_lo ) then
2466  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2467  else
2468  ftind = zr1 / (bpade + zr1)
2469  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2470  zexp4 = exp_tbl(itind)
2471  endif
2472 
2473  zldbt0(k) = zexp4
2474  ztdbt0 = zexp4 * ztdbt0
2475  enddo ! end do_k_loop
2476 
2477  call swflux &
2478 ! --- inputs:
2479  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2480  & nlay, nlp1, &
2481 ! --- outputs:
2482  & zfu, zfd &
2483  & )
2484 
2485 ! --- ... compute upward and downward fluxes at levels
2486  do k = 1, nlp1
2487  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
2488  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
2489  enddo
2490 
2491 !! --- ... surface downward beam/diffused flux components
2492  zb1 = zsolar*ztdbt0
2493  zb2 = zsolar*(zfd(1) - ztdbt0)
2494 
2495  if (ibd /= 0) then
2496  sfbm0(ibd) = sfbm0(ibd) + zb1
2497  sfdf0(ibd) = sfdf0(ibd) + zb2
2498  else
2499  zf1 = 0.5 * zb1
2500  zf2 = 0.5 * zb2
2501  sfbm0(1) = sfbm0(1) + zf1
2502  sfdf0(1) = sfdf0(1) + zf2
2503  sfbm0(2) = sfbm0(2) + zf1
2504  sfdf0(2) = sfdf0(2) + zf2
2505  endif
2506 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
2507 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
2508 
2509 ! --- ... compute total sky optical parameters, layer reflectance and transmittance
2510 
2511  if ( cf1 > eps ) then
2512 
2513 ! --- ... set up toa direct beam and surface values (beam and diff)
2514  ztdbt0 = f_one
2515  zldbt(1) = f_zero
2516 
2517  do k = nlay, 1, -1
2518  kp = k + 1
2519  zc0 = f_one - cldfrc(k)
2520  zc1 = cldfrc(k)
2521  if ( zc1 > ftiny ) then ! it is a cloudy-layer
2522 
2523  ztau0 = ztaus(k) + taucw(k,ib)
2524  zssa0 = zssas(k) + ssacw(k,ib)
2525  zasy0 = zasys(k) + asycw(k,ib)
2526  zssaw = min(oneminus, zssa0 / ztau0)
2527  zasyw = zasy0 / max(ftiny, zssa0)
2528 
2529 ! --- ... delta scaling for total-sky condition
2530  za1 = zasyw * zasyw
2531  za2 = zssaw * za1
2532 
2533  ztau1 = (f_one - za2) * ztau0
2534  zssa1 = (zssaw - za2) / (f_one - za2)
2535 !org zasy1 = (zasyw - za1) / (f_one - za1)
2536  zasy1 = zasyw / (f_one + zasyw)
2537  zasy3 = 0.75 * zasy1
2538 
2539 ! --- ... general two-stream expressions
2540  if ( iswmode == 1 ) then
2541  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
2542  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
2543  zgam3 = 0.5 - zasy3 * cosz
2544  elseif ( iswmode == 2 ) then ! pifm
2545  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
2546  zgam2 = 0.75* zssa1 * (f_one- zasy1)
2547  zgam3 = 0.5 - zasy3 * cosz
2548  elseif ( iswmode == 3 ) then ! discrete ordinates
2549  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
2550  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
2551  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
2552  endif
2553  zgam4 = f_one - zgam3
2554 
2555  zrefb1 = zrefb(kp)
2556  zrefd1 = zrefd(kp)
2557  ztrab1 = ztrab(kp)
2558  ztrad1 = ztrad(kp)
2559 
2560 ! --- ... compute homogeneous reflectance and transmittance
2561 
2562  if ( zssaw >= zcrit ) then ! for conservative scattering
2563  za1 = zgam1 * cosz - zgam3
2564  za2 = zgam1 * ztau1
2565 
2566 ! --- ... use exponential lookup table for transmittance, or expansion
2567 ! of exponential for low optical depth
2568 
2569  zb1 = min( ztau1*sntz , 500.0 )
2570  if ( zb1 <= od_lo ) then
2571  zb2 = f_one - zb1 + 0.5*zb1*zb1
2572  else
2573  ftind = zb1 / (bpade + zb1)
2574  itind = ftind*ntbmx + 0.5
2575  zb2 = exp_tbl(itind)
2576  endif
2577 
2578 ! ... collimated beam
2579  zrefb(kp) = max(f_zero, min(f_one, &
2580  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
2581  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
2582 
2583 ! ... isotropic incidence
2584  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
2585  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
2586 
2587  else ! for non-conservative scattering
2588  za1 = zgam1*zgam4 + zgam2*zgam3
2589  za2 = zgam1*zgam3 + zgam2*zgam4
2590  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
2591  zrk2= 2.0 * zrk
2592 
2593  zrp = zrk * cosz
2594  zrp1 = f_one + zrp
2595  zrm1 = f_one - zrp
2596  zrpp = f_one - zrp*zrp
2597  zrkg1= zrk + zgam1
2598  zrkg3= zrk * zgam3
2599  zrkg4= zrk * zgam4
2600 
2601  zr1 = zrm1 * (za2 + zrkg3)
2602  zr2 = zrp1 * (za2 - zrkg3)
2603  zr3 = zrk2 * (zgam3 - za2*cosz)
2604  zr4 = zrpp * zrkg1
2605  zr5 = zrpp * (zrk - zgam1)
2606 
2607  zt1 = zrp1 * (za1 + zrkg4)
2608  zt2 = zrm1 * (za1 - zrkg4)
2609  zt3 = zrk2 * (zgam4 + za1*cosz)
2610 
2611 ! --- ... use exponential lookup table for transmittance, or expansion
2612 ! of exponential for low optical depth
2613 
2614  zb1 = min( zrk*ztau1, 500.0 )
2615  if ( zb1 <= od_lo ) then
2616  zexm1 = f_one - zb1 + 0.5*zb1*zb1
2617  else
2618  ftind = zb1 / (bpade + zb1)
2619  itind = ftind*ntbmx + 0.5
2620  zexm1 = exp_tbl(itind)
2621  endif
2622  zexp1 = f_one / zexm1
2623 
2624  zb2 = min( ztau1*sntz, 500.0 )
2625  if ( zb2 <= od_lo ) then
2626  zexm2 = f_one - zb2 + 0.5*zb2*zb2
2627  else
2628  ftind = zb2 / (bpade + zb2)
2629  itind = ftind*ntbmx + 0.5
2630  zexm2 = exp_tbl(itind)
2631  endif
2632  zexp2 = f_one / zexm2
2633  ze1r45 = zr4*zexp1 + zr5*zexm1
2634 
2635 ! ... collimated beam
2636  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
2637  zrefb(kp) = eps1
2638  ztrab(kp) = zexm2
2639  else
2640  zden1 = zssa1 / ze1r45
2641  zrefb(kp) = max(f_zero, min(f_one, &
2642  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
2643  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
2644  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
2645  endif
2646 
2647 ! ... diffuse beam
2648  zden1 = zr4 / (ze1r45 * zrkg1)
2649  zrefd(kp) = max(f_zero, min(f_one, &
2650  & zgam2*(zexp1 - zexm1)*zden1 ))
2651  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
2652  endif ! end if_zssaw_block
2653 
2654 ! --- ... combine clear and cloudy contributions for total sky
2655 ! and calculate direct beam transmittances
2656 
2657  zrefb(kp) = zc0*zrefb1 + zc1*zrefb(kp)
2658  zrefd(kp) = zc0*zrefd1 + zc1*zrefd(kp)
2659  ztrab(kp) = zc0*ztrab1 + zc1*ztrab(kp)
2660  ztrad(kp) = zc0*ztrad1 + zc1*ztrad(kp)
2661 
2662 ! --- ... direct beam transmittance. use exponential lookup table
2663 ! for transmittance, or expansion of exponential for low
2664 ! optical depth
2665 
2666  zr1 = ztau1 * sntz
2667  if ( zr1 <= od_lo ) then
2668  zexp3 = f_one - zr1 + 0.5*zr1*zr1
2669  else
2670  ftind = zr1 / (bpade + zr1)
2671  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2672  zexp3 = exp_tbl(itind)
2673  endif
2674 
2675  zldbt(kp) = zc0*zldbt(kp) + zc1*zexp3
2676  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2677 
2678 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2679 ! (must use 'orig', unscaled cloud optical depth)
2680 
2681  zr1 = ztau0 * sntz
2682  if ( zr1 <= od_lo ) then
2683  zexp4 = f_one - zr1 + 0.5*zr1*zr1
2684  else
2685  ftind = zr1 / (bpade + zr1)
2686  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
2687  zexp4 = exp_tbl(itind)
2688  endif
2689 
2690  ztdbt0 = (zc0*zldbt0(k) + zc1*zexp4) * ztdbt0
2691 
2692  else ! if_zc1_block --- it is a clear layer
2693 
2694 ! --- ... direct beam transmittance
2695  ztdbt(k) = zldbt(kp) * ztdbt(kp)
2696 
2697 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
2698  ztdbt0 = zldbt0(k) * ztdbt0
2699 
2700  endif ! end if_zc1_block
2701  enddo ! end do_k_loop
2702 
2703 ! --- ... perform vertical quadrature
2704 
2705  call swflux &
2706 ! --- inputs:
2707  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
2708  & nlay, nlp1, &
2709 ! --- outputs:
2710  & zfu, zfd &
2711  & )
2712 
2713 ! --- ... compute upward and downward fluxes at levels
2714  do k = 1, nlp1
2715  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
2716  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
2717  enddo
2718 
2719 !! --- ... surface downward beam/diffused flux components
2720  zb1 = zsolar*ztdbt0
2721  zb2 = zsolar*(zfd(1) - ztdbt0)
2722 
2723  if (ibd /= 0) then
2724  sfbmc(ibd) = sfbmc(ibd) + zb1
2725  sfdfc(ibd) = sfdfc(ibd) + zb2
2726  else
2727  zf1 = 0.5 * zb1
2728  zf2 = 0.5 * zb2
2729  sfbmc(1) = sfbmc(1) + zf1
2730  sfdfc(1) = sfdfc(1) + zf2
2731  sfbmc(2) = sfbmc(2) + zf1
2732  sfdfc(2) = sfdfc(2) + zf2
2733  endif
2734 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
2735 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
2736 
2737  endif ! end if_cf1_block
2738 
2739  enddo lab_do_jg
2740 
2741 ! --- ... end of g-point loop
2742 
2743  do ib = 1, nbdsw
2744  ftoadc = ftoadc + fxdn0(nlp1,ib)
2745  ftoau0 = ftoau0 + fxup0(nlp1,ib)
2746  fsfcu0 = fsfcu0 + fxup0(1,ib)
2747  fsfcd0 = fsfcd0 + fxdn0(1,ib)
2748  enddo
2749 
2750 !! --- ... uv-b surface downward flux
2751  ibd = nuvb - nblow + 1
2752  suvbf0 = fxdn0(1,ibd)
2753 
2754  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
2755  do ib = 1, nbdsw
2756  do k = 1, nlp1
2757  fxupc(k,ib) = fxup0(k,ib)
2758  fxdnc(k,ib) = fxdn0(k,ib)
2759  enddo
2760  enddo
2761 
2762  ftoauc = ftoau0
2763  fsfcuc = fsfcu0
2764  fsfcdc = fsfcd0
2765 
2766 !! --- ... surface downward beam/diffused flux components
2767  sfbmc(1) = sfbm0(1)
2768  sfdfc(1) = sfdf0(1)
2769  sfbmc(2) = sfbm0(2)
2770  sfdfc(2) = sfdf0(2)
2771 
2772 !! --- ... uv-b surface downward flux
2773  suvbfc = suvbf0
2774  else ! cloudy column, compute total-sky fluxes
2775  do ib = 1, nbdsw
2776  do k = 1, nlp1
2777  fxupc(k,ib) = cf1*fxupc(k,ib) + cf0*fxup0(k,ib)
2778  fxdnc(k,ib) = cf1*fxdnc(k,ib) + cf0*fxdn0(k,ib)
2779  enddo
2780  enddo
2781 
2782  do ib = 1, nbdsw
2783  ftoauc = ftoauc + fxupc(nlp1,ib)
2784  fsfcuc = fsfcuc + fxupc(1,ib)
2785  fsfcdc = fsfcdc + fxdnc(1,ib)
2786  enddo
2787 
2788 !! --- ... uv-b surface downward flux
2789  suvbfc = fxdnc(1,ibd)
2790 
2791 !! --- ... surface downward beam/diffused flux components
2792  sfbmc(1) = cf1*sfbmc(1) + cf0*sfbm0(1)
2793  sfbmc(2) = cf1*sfbmc(2) + cf0*sfbm0(2)
2794  sfdfc(1) = cf1*sfdfc(1) + cf0*sfdf0(1)
2795  sfdfc(2) = cf1*sfdfc(2) + cf0*sfdf0(2)
2796  endif ! end if_cf1_block
2797 
2798  return
2799 !...................................

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine module_radsw_main::spcvrtm ( )
private

Definition at line 2809 of file radsw_main.f.

References bpade, eps, exp_tbl, f_one, f_zero, ftiny, idxsfc, physparam::iswmode, module_radsw_parameters::nbdsw, module_radsw_parameters::nblow, module_radsw_parameters::ngb, module_radsw_parameters::ngptsw, module_radsw_parameters::ntbmx, nuvb, oneminus, and swflux().

Referenced by swrad().

2809 ! --- inputs:
2810  & ( ssolar,cosz,sntz,albbm,albdf,sfluxzen,cldfmc, &
2811  & cf1,cf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
2812  & nlay, nlp1, &
2813 ! --- outputs:
2814  & fxupc,fxdnc,fxup0,fxdn0, &
2815  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
2816  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
2817  & )
2818 
2819 ! =================== program usage description =================== !
2820 ! !
2821 ! purpose: computes the shortwave radiative fluxes using two-stream !
2822 ! method of h. barker and mcica, the monte-carlo independent!
2823 ! column approximation, for the representation of sub-grid !
2824 ! cloud variability (i.e. cloud overlap). !
2825 ! !
2826 ! subprograms called: swflux !
2827 ! !
2828 ! ==================== defination of variables ==================== !
2829 ! !
2830 ! inputs: size !
2831 ! ssolar - real, incoming solar flux at top 1 !
2832 ! cosz - real, cosine solar zenith angle 1 !
2833 ! sntz - real, secant solar zenith angle 1 !
2834 ! albbm - real, surface albedo for direct beam radiation 2 !
2835 ! albdf - real, surface albedo for diffused radiation 2 !
2836 ! sfluxzen- real, spectral distribution of incoming solar flux ngptsw!
2837 ! cldfmc - real, layer cloud fraction for g-point nlay*ngptsw!
2838 ! cf1 - real, >0: cloudy sky, otherwise: clear sky 1 !
2839 ! cf0 - real, =1-cf1 1 !
2840 ! taug - real, spectral optical depth for gases nlay*ngptsw!
2841 ! taur - real, optical depth for rayleigh scattering nlay*ngptsw!
2842 ! tauae - real, aerosols optical depth nlay*nbdsw !
2843 ! ssaae - real, aerosols single scattering albedo nlay*nbdsw !
2844 ! asyae - real, aerosols asymmetry factor nlay*nbdsw !
2845 ! taucw - real, weighted cloud optical depth nlay*nbdsw !
2846 ! ssacw - real, weighted cloud single scat albedo nlay*nbdsw !
2847 ! asycw - real, weighted cloud asymmetry factor nlay*nbdsw !
2848 ! nlay,nlp1 - integer, number of layers/levels 1 !
2849 ! !
2850 ! output variables: !
2851 ! fxupc - real, tot sky upward flux nlp1*nbdsw !
2852 ! fxdnc - real, tot sky downward flux nlp1*nbdsw !
2853 ! fxup0 - real, clr sky upward flux nlp1*nbdsw !
2854 ! fxdn0 - real, clr sky downward flux nlp1*nbdsw !
2855 ! ftoauc - real, tot sky toa upwd flux 1 !
2856 ! ftoau0 - real, clr sky toa upwd flux 1 !
2857 ! ftoadc - real, toa downward (incoming) solar flux 1 !
2858 ! fsfcuc - real, tot sky sfc upwd flux 1 !
2859 ! fsfcu0 - real, clr sky sfc upwd flux 1 !
2860 ! fsfcdc - real, tot sky sfc dnwd flux 1 !
2861 ! fsfcd0 - real, clr sky sfc dnwd flux 1 !
2862 ! sfbmc - real, tot sky sfc dnwd beam flux (nir/uv+vis) 2 !
2863 ! sfdfc - real, tot sky sfc dnwd diff flux (nir/uv+vis) 2 !
2864 ! sfbm0 - real, clr sky sfc dnwd beam flux (nir/uv+vis) 2 !
2865 ! sfdf0 - real, clr sky sfc dnwd diff flux (nir/uv+vis) 2 !
2866 ! suvbfc - real, tot sky sfc dnwd uv-b flux 1 !
2867 ! suvbf0 - real, clr sky sfc dnwd uv-b flux 1 !
2868 ! !
2869 ! internal variables: !
2870 ! zrefb - real, direct beam reflectivity for clear/cloudy nlp1 !
2871 ! zrefd - real, diffuse reflectivity for clear/cloudy nlp1 !
2872 ! ztrab - real, direct beam transmissivity for clear/cloudy nlp1 !
2873 ! ztrad - real, diffuse transmissivity for clear/cloudy nlp1 !
2874 ! zldbt - real, layer beam transmittance for clear/cloudy nlp1 !
2875 ! ztdbt - real, lev total beam transmittance for clr/cld nlp1 !
2876 ! !
2877 ! control parameters in module "physparam" !
2878 ! iswmode - control flag for 2-stream transfer schemes !
2879 ! = 1 delta-eddington (joseph et al., 1976) !
2880 ! = 2 pifm (zdunkowski et al., 1980) !
2881 ! = 3 discrete ordinates (liou, 1973) !
2882 ! !
2883 ! ******************************************************************* !
2884 ! original code description !
2885 ! !
2886 ! method: !
2887 ! ------- !
2888 ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations. !
2889 ! kmodts = 1 eddington (joseph et al., 1976) !
2890 ! = 2 pifm (zdunkowski et al., 1980) !
2891 ! = 3 discrete ordinates (liou, 1973) !
2892 ! !
2893 ! modifications: !
2894 ! -------------- !
2895 ! original: h. barker !
2896 ! revision: merge with rrtmg_sw: j.-j.morcrette, ecmwf, feb 2003 !
2897 ! revision: add adjustment for earth/sun distance:mjiacono,aer,oct2003!
2898 ! revision: bug fix for use of palbp and palbd: mjiacono, aer, nov2003!
2899 ! revision: bug fix to apply delta scaling to clear sky: aer, dec2004 !
2900 ! revision: code modified so that delta scaling is not done in cloudy !
2901 ! profiles if routine cldprop is used; delta scaling can be !
2902 ! applied by swithcing code below if cldprop is not used to !
2903 ! get cloud properties. aer, jan 2005 !
2904 ! revision: uniform formatting for rrtmg: mjiacono, aer, jul 2006 !
2905 ! revision: use exponential lookup table for transmittance: mjiacono, !
2906 ! aer, aug 2007 !
2907 ! !
2908 ! ******************************************************************* !
2909 ! ====================== end of description block ================= !
2910 
2911 ! --- constant parameters:
2912  real (kind=kind_phys), parameter :: zcrit = 0.9999995 ! thresold for conservative scattering
2913  real (kind=kind_phys), parameter :: zsr3 = sqrt(3.0)
2914  real (kind=kind_phys), parameter :: od_lo = 0.06
2915  real (kind=kind_phys), parameter :: eps1 = 1.0e-8
2916 
2917 ! --- inputs:
2918  integer, intent(in) :: nlay, nlp1
2919 
2920  real (kind=kind_phys), dimension(nlay,ngptsw), intent(in) :: &
2921  & taug, taur, cldfmc
2922  real (kind=kind_phys), dimension(nlay,nbdsw), intent(in) :: &
2923  & taucw, ssacw, asycw, tauae, ssaae, asyae
2924 
2925  real (kind=kind_phys), dimension(ngptsw), intent(in) :: sfluxzen
2926 
2927  real (kind=kind_phys), dimension(2), intent(in) :: albbm, albdf
2928 
2929  real (kind=kind_phys), intent(in) :: cosz, sntz, cf1, cf0, ssolar
2930 
2931 ! --- outputs:
2932  real (kind=kind_phys), dimension(nlp1,nbdsw), intent(out) :: &
2933  & fxupc, fxdnc, fxup0, fxdn0
2934 
2935  real (kind=kind_phys), dimension(2), intent(out) :: sfbmc, sfdfc, &
2936  & sfbm0, sfdf0
2937 
2938  real (kind=kind_phys), intent(out) :: suvbfc, suvbf0, ftoadc, &
2939  & ftoauc, ftoau0, fsfcuc, fsfcu0, fsfcdc, fsfcd0
2940 
2941 ! --- locals:
2942  real (kind=kind_phys), dimension(nlay) :: ztaus, zssas, zasys, &
2943  & zldbt0
2944 
2945  real (kind=kind_phys), dimension(nlp1) :: zrefb, zrefd, ztrab, &
2946  & ztrad, ztdbt, zldbt, zfu, zfd
2947 
2948  real (kind=kind_phys) :: ztau1, zssa1, zasy1, ztau0, zssa0, &
2949  & zasy0, zasy3, zssaw, zasyw, zgam1, zgam2, zgam3, zgam4, &
2950  & za1, za2, zb1, zb2, zrk, zrk2, zrp, zrp1, zrm1, zrpp, &
2951  & zrkg1, zrkg3, zrkg4, zexp1, zexm1, zexp2, zexm2, zden1, &
2952  & zexp3, zexp4, ze1r45, ftind, zsolar, ztdbt0, zr1, zr2, &
2953  & zr3, zr4, zr5, zt1, zt2, zt3, zf1, zf2
2954 
2955  integer :: ib, ibd, jb, jg, k, kp, itind
2956 !
2957 !===> ... begin here
2958 !
2959 ! --- ... initialization of output fluxes
2960 
2961  do ib = 1, nbdsw
2962  do k = 1, nlp1
2963  fxdnc(k,ib) = f_zero
2964  fxupc(k,ib) = f_zero
2965  fxdn0(k,ib) = f_zero
2966  fxup0(k,ib) = f_zero
2967  enddo
2968  enddo
2969 
2970  ftoadc = f_zero
2971  ftoauc = f_zero
2972  ftoau0 = f_zero
2973  fsfcuc = f_zero
2974  fsfcu0 = f_zero
2975  fsfcdc = f_zero
2976  fsfcd0 = f_zero
2977 
2978 !! --- ... uv-b surface downward fluxes
2979  suvbfc = f_zero
2980  suvbf0 = f_zero
2981 
2982 !! --- ... output surface flux components
2983  sfbmc(1) = f_zero
2984  sfbmc(2) = f_zero
2985  sfdfc(1) = f_zero
2986  sfdfc(2) = f_zero
2987  sfbm0(1) = f_zero
2988  sfbm0(2) = f_zero
2989  sfdf0(1) = f_zero
2990  sfdf0(2) = f_zero
2991 
2992 ! --- ... loop over all g-points in each band
2993 
2994  lab_do_jg : do jg = 1, ngptsw
2995 
2996  jb = ngb(jg)
2997  ib = jb + 1 - nblow
2998  ibd = idxsfc(jb) ! spectral band index
2999 
3000  zsolar = ssolar * sfluxzen(jg)
3001 
3002 ! --- ... set up toa direct beam and surface values (beam and diff)
3003 
3004  ztdbt(nlp1) = f_one
3005  ztdbt0 = f_one
3006 
3007  zldbt(1) = f_zero
3008  if (ibd /= 0) then
3009  zrefb(1) = albbm(ibd)
3010  zrefd(1) = albdf(ibd)
3011  else
3012  zrefb(1) = 0.5 * (albbm(1) + albbm(2))
3013  zrefd(1) = 0.5 * (albdf(1) + albdf(2))
3014  endif
3015  ztrab(1) = f_zero
3016  ztrad(1) = f_zero
3017 
3018 ! --- ... compute clear-sky optical parameters, layer reflectance and transmittance
3019 
3020  do k = nlay, 1, -1
3021  kp = k + 1
3022 
3023  ztau0 = max( ftiny, taur(k,jg)+taug(k,jg)+tauae(k,ib) )
3024  zssa0 = taur(k,jg) + tauae(k,ib)*ssaae(k,ib)
3025  zasy0 = asyae(k,ib)*ssaae(k,ib)*tauae(k,ib)
3026  zssaw = min( oneminus, zssa0 / ztau0 )
3027  zasyw = zasy0 / max( ftiny, zssa0 )
3028 
3029 ! --- ... saving clear-sky quantities for later total-sky usage
3030  ztaus(k) = ztau0
3031  zssas(k) = zssa0
3032  zasys(k) = zasy0
3033 
3034 ! --- ... delta scaling for clear-sky condition
3035  za1 = zasyw * zasyw
3036  za2 = zssaw * za1
3037 
3038  ztau1 = (f_one - za2) * ztau0
3039  zssa1 = (zssaw - za2) / (f_one - za2)
3040 !org zasy1 = (zasyw - za1) / (f_one - za1) ! this line is replaced by the next
3041  zasy1 = zasyw / (f_one + zasyw) ! to reduce truncation error
3042  zasy3 = 0.75 * zasy1
3043 
3044 ! --- ... general two-stream expressions
3045  if ( iswmode == 1 ) then
3046  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3047  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
3048  zgam3 = 0.5 - zasy3 * cosz
3049  elseif ( iswmode == 2 ) then ! pifm
3050  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3051  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3052  zgam3 = 0.5 - zasy3 * cosz
3053  elseif ( iswmode == 3 ) then ! discrete ordinates
3054  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3055  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3056  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3057  endif
3058  zgam4 = f_one - zgam3
3059 
3060 ! --- ... compute homogeneous reflectance and transmittance
3061 
3062  if ( zssaw >= zcrit ) then ! for conservative scattering
3063  za1 = zgam1 * cosz - zgam3
3064  za2 = zgam1 * ztau1
3065 
3066 ! --- ... use exponential lookup table for transmittance, or expansion
3067 ! of exponential for low optical depth
3068 
3069  zb1 = min( ztau1*sntz , 500.0 )
3070  if ( zb1 <= od_lo ) then
3071  zb2 = f_one - zb1 + 0.5*zb1*zb1
3072  else
3073  ftind = zb1 / (bpade + zb1)
3074  itind = ftind*ntbmx + 0.5
3075  zb2 = exp_tbl(itind)
3076  endif
3077 
3078 ! ... collimated beam
3079  zrefb(kp) = max(f_zero, min(f_one, &
3080  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3081  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp) ))
3082 
3083 ! ... isotropic incidence
3084  zrefd(kp) = max(f_zero, min(f_one, za2/(f_one + za2) ))
3085  ztrad(kp) = max(f_zero, min(f_one, f_one-zrefd(kp) ))
3086 
3087  else ! for non-conservative scattering
3088  za1 = zgam1*zgam4 + zgam2*zgam3
3089  za2 = zgam1*zgam3 + zgam2*zgam4
3090  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3091  zrk2= 2.0 * zrk
3092 
3093  zrp = zrk * cosz
3094  zrp1 = f_one + zrp
3095  zrm1 = f_one - zrp
3096  zrpp = f_one - zrp*zrp
3097  zrkg1= zrk + zgam1
3098  zrkg3= zrk * zgam3
3099  zrkg4= zrk * zgam4
3100 
3101  zr1 = zrm1 * (za2 + zrkg3)
3102  zr2 = zrp1 * (za2 - zrkg3)
3103  zr3 = zrk2 * (zgam3 - za2*cosz)
3104  zr4 = zrpp * zrkg1
3105  zr5 = zrpp * (zrk - zgam1)
3106 
3107  zt1 = zrp1 * (za1 + zrkg4)
3108  zt2 = zrm1 * (za1 - zrkg4)
3109  zt3 = zrk2 * (zgam4 + za1*cosz)
3110 
3111 ! --- ... use exponential lookup table for transmittance, or expansion
3112 ! of exponential for low optical depth
3113 
3114  zb1 = min( zrk*ztau1, 500.0 )
3115  if ( zb1 <= od_lo ) then
3116  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3117  else
3118  ftind = zb1 / (bpade + zb1)
3119  itind = ftind*ntbmx + 0.5
3120  zexm1 = exp_tbl(itind)
3121  endif
3122  zexp1 = f_one / zexm1
3123 
3124  zb2 = min( sntz*ztau1, 500.0 )
3125  if ( zb2 <= od_lo ) then
3126  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3127  else
3128  ftind = zb2 / (bpade + zb2)
3129  itind = ftind*ntbmx + 0.5
3130  zexm2 = exp_tbl(itind)
3131  endif
3132  zexp2 = f_one / zexm2
3133  ze1r45 = zr4*zexp1 + zr5*zexm1
3134 
3135 ! ... collimated beam
3136  if (ze1r45>=-eps1 .and. ze1r45<=eps1) then
3137  zrefb(kp) = eps1
3138  ztrab(kp) = zexm2
3139  else
3140  zden1 = zssa1 / ze1r45
3141  zrefb(kp) = max(f_zero, min(f_one, &
3142  & (zr1*zexp1 - zr2*zexm1 - zr3*zexm2)*zden1 ))
3143  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one &
3144  & - (zt1*zexp1 - zt2*zexm1 - zt3*zexp2)*zden1) ))
3145  endif
3146 
3147 ! ... diffuse beam
3148  zden1 = zr4 / (ze1r45 * zrkg1)
3149  zrefd(kp) = max(f_zero, min(f_one, &
3150  & zgam2*(zexp1 - zexm1)*zden1 ))
3151  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3152  endif ! end if_zssaw_block
3153 
3154 ! --- ... direct beam transmittance. use exponential lookup table
3155 ! for transmittance, or expansion of exponential for low
3156 ! optical depth
3157 
3158  zr1 = ztau1 * sntz
3159  if ( zr1 <= od_lo ) then
3160  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3161  else
3162  ftind = zr1 / (bpade + zr1)
3163  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3164  zexp3 = exp_tbl(itind)
3165  endif
3166 
3167  ztdbt(k) = zexp3 * ztdbt(kp)
3168  zldbt(kp) = zexp3
3169 
3170 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3171 ! (must use 'orig', unscaled cloud optical depth)
3172 
3173  zr1 = ztau0 * sntz
3174  if ( zr1 <= od_lo ) then
3175  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3176  else
3177  ftind = zr1 / (bpade + zr1)
3178  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3179  zexp4 = exp_tbl(itind)
3180  endif
3181 
3182  zldbt0(k) = zexp4
3183  ztdbt0 = zexp4 * ztdbt0
3184  enddo ! end do_k_loop
3185 
3186  call swflux &
3187 ! --- inputs:
3188  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3189  & nlay, nlp1, &
3190 ! --- outputs:
3191  & zfu, zfd &
3192  & )
3193 
3194 ! --- ... compute upward and downward fluxes at levels
3195  do k = 1, nlp1
3196  fxup0(k,ib) = fxup0(k,ib) + zsolar*zfu(k)
3197  fxdn0(k,ib) = fxdn0(k,ib) + zsolar*zfd(k)
3198  enddo
3199 
3200 !! --- ... surface downward beam/diffuse flux components
3201  zb1 = zsolar*ztdbt0
3202  zb2 = zsolar*(zfd(1) - ztdbt0)
3203 
3204  if (ibd /= 0) then
3205  sfbm0(ibd) = sfbm0(ibd) + zb1
3206  sfdf0(ibd) = sfdf0(ibd) + zb2
3207  else
3208  zf1 = 0.5 * zb1
3209  zf2 = 0.5 * zb2
3210  sfbm0(1) = sfbm0(1) + zf1
3211  sfdf0(1) = sfdf0(1) + zf2
3212  sfbm0(2) = sfbm0(2) + zf1
3213  sfdf0(2) = sfdf0(2) + zf2
3214  endif
3215 ! sfbm0(ibd) = sfbm0(ibd) + zsolar*ztdbt0
3216 ! sfdf0(ibd) = sfdf0(ibd) + zsolar*(zfd(1) - ztdbt0)
3217 
3218 ! --- ... compute total sky optical parameters, layer reflectance and transmittance
3219 
3220  if ( cf1 > eps ) then
3221 
3222 ! --- ... set up toa direct beam and surface values (beam and diff)
3223  ztdbt0 = f_one
3224  zldbt(1) = f_zero
3225 
3226  do k = nlay, 1, -1
3227  kp = k + 1
3228  if ( cldfmc(k,jg) > ftiny ) then ! it is a cloudy-layer
3229 
3230  ztau0 = ztaus(k) + taucw(k,ib)
3231  zssa0 = zssas(k) + ssacw(k,ib)
3232  zasy0 = zasys(k) + asycw(k,ib)
3233  zssaw = min(oneminus, zssa0 / ztau0)
3234  zasyw = zasy0 / max(ftiny, zssa0)
3235 
3236 ! --- ... delta scaling for total-sky condition
3237  za1 = zasyw * zasyw
3238  za2 = zssaw * za1
3239 
3240  ztau1 = (f_one - za2) * ztau0
3241  zssa1 = (zssaw - za2) / (f_one - za2)
3242 !org zasy1 = (zasyw - za1) / (f_one - za1)
3243  zasy1 = zasyw / (f_one + zasyw)
3244  zasy3 = 0.75 * zasy1
3245 
3246 ! --- ... general two-stream expressions
3247  if ( iswmode == 1 ) then
3248  zgam1 = 1.75 - zssa1 * (f_one + zasy3)
3249  zgam2 =-0.25 - zssa1 * (f_one - zasy3)
3250  zgam3 = 0.5 - zasy3 * cosz
3251  elseif ( iswmode == 2 ) then ! pifm
3252  zgam1 = 2.0 - zssa1 * (1.25 + zasy3)
3253  zgam2 = 0.75* zssa1 * (f_one- zasy1)
3254  zgam3 = 0.5 - zasy3 * cosz
3255  elseif ( iswmode == 3 ) then ! discrete ordinates
3256  zgam1 = zsr3 * (2.0 - zssa1 * (1.0 + zasy1)) * 0.5
3257  zgam2 = zsr3 * zssa1 * (1.0 - zasy1) * 0.5
3258  zgam3 = (1.0 - zsr3 * zasy1 * cosz) * 0.5
3259  endif
3260  zgam4 = f_one - zgam3
3261 
3262 ! --- ... compute homogeneous reflectance and transmittance
3263 
3264  if ( zssaw >= zcrit ) then ! for conservative scattering
3265  za1 = zgam1 * cosz - zgam3
3266  za2 = zgam1 * ztau1
3267 
3268 ! --- ... use exponential lookup table for transmittance, or expansion
3269 ! of exponential for low optical depth
3270 
3271  zb1 = min( ztau1*sntz , 500.0 )
3272  if ( zb1 <= od_lo ) then
3273  zb2 = f_one - zb1 + 0.5*zb1*zb1
3274  else
3275  ftind = zb1 / (bpade + zb1)
3276  itind = ftind*ntbmx + 0.5
3277  zb2 = exp_tbl(itind)
3278  endif
3279 
3280 ! ... collimated beam
3281  zrefb(kp) = max(f_zero, min(f_one, &
3282  & (za2 - za1*(f_one - zb2))/(f_one + za2) ))
3283  ztrab(kp) = max(f_zero, min(f_one, f_one-zrefb(kp)))
3284 
3285 ! ... isotropic incidence
3286  zrefd(kp) = max(f_zero, min(f_one, za2 / (f_one+za2) ))
3287  ztrad(kp) = max(f_zero, min(f_one, f_one - zrefd(kp) ))
3288 
3289  else ! for non-conservative scattering
3290  za1 = zgam1*zgam4 + zgam2*zgam3
3291  za2 = zgam1*zgam3 + zgam2*zgam4
3292  zrk = sqrt( (zgam1 - zgam2) * (zgam1 + zgam2) )
3293  zrk2= 2.0 * zrk
3294 
3295  zrp = zrk * cosz
3296  zrp1 = f_one + zrp
3297  zrm1 = f_one - zrp
3298  zrpp = f_one - zrp*zrp
3299  zrkg1= zrk + zgam1
3300  zrkg3= zrk * zgam3
3301  zrkg4= zrk * zgam4
3302 
3303  zr1 = zrm1 * (za2 + zrkg3)
3304  zr2 = zrp1 * (za2 - zrkg3)
3305  zr3 = zrk2 * (zgam3 - za2*cosz)
3306  zr4 = zrpp * zrkg1
3307  zr5 = zrpp * (zrk - zgam1)
3308 
3309  zt1 = zrp1 * (za1 + zrkg4)
3310  zt2 = zrm1 * (za1 - zrkg4)
3311  zt3 = zrk2 * (zgam4 + za1*cosz)
3312 
3313 ! --- ... use exponential lookup table for transmittance, or expansion
3314 ! of exponential for low optical depth
3315 
3316  zb1 = min( zrk*ztau1, 500.0 )
3317  if ( zb1 <= od_lo ) then
3318  zexm1 = f_one - zb1 + 0.5*zb1*zb1
3319  else
3320  ftind = zb1 / (bpade + zb1)
3321  itind = ftind*ntbmx + 0.5
3322  zexm1 = exp_tbl(itind)
3323  endif
3324  zexp1 = f_one / zexm1
3325 
3326  zb2 = min( ztau1*sntz, 500.0 )
3327  if ( zb2 <= od_lo ) then
3328  zexm2 = f_one - zb2 + 0.5*zb2*zb2
3329  else
3330  ftind = zb2 / (bpade + zb2)
3331  itind = ftind*ntbmx + 0.5
3332  zexm2 = exp_tbl(itind)
3333  endif
3334  zexp2 = f_one / zexm2
3335  ze1r45 = zr4*zexp1 + zr5*zexm1
3336 
3337 ! ... collimated beam
3338  if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then
3339  zrefb(kp) = eps1
3340  ztrab(kp) = zexm2
3341  else
3342  zden1 = zssa1 / ze1r45
3343  zrefb(kp) = max(f_zero, min(f_one, &
3344  & (zr1*zexp1-zr2*zexm1-zr3*zexm2)*zden1 ))
3345  ztrab(kp) = max(f_zero, min(f_one, zexm2*(f_one - &
3346  & (zt1*zexp1-zt2*zexm1-zt3*zexp2)*zden1) ))
3347  endif
3348 
3349 ! ... diffuse beam
3350  zden1 = zr4 / (ze1r45 * zrkg1)
3351  zrefd(kp) = max(f_zero, min(f_one, &
3352  & zgam2*(zexp1 - zexm1)*zden1 ))
3353  ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 ))
3354  endif ! end if_zssaw_block
3355 
3356 ! --- ... direct beam transmittance. use exponential lookup table
3357 ! for transmittance, or expansion of exponential for low
3358 ! optical depth
3359 
3360  zr1 = ztau1 * sntz
3361  if ( zr1 <= od_lo ) then
3362  zexp3 = f_one - zr1 + 0.5*zr1*zr1
3363  else
3364  ftind = zr1 / (bpade + zr1)
3365  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3366  zexp3 = exp_tbl(itind)
3367  endif
3368 
3369  zldbt(kp) = zexp3
3370  ztdbt(k) = zexp3 * ztdbt(kp)
3371 
3372 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3373 ! (must use 'orig', unscaled cloud optical depth)
3374 
3375  zr1 = ztau0 * sntz
3376  if ( zr1 <= od_lo ) then
3377  zexp4 = f_one - zr1 + 0.5*zr1*zr1
3378  else
3379  ftind = zr1 / (bpade + zr1)
3380  itind = max(0, min(ntbmx, int(0.5+ntbmx*ftind) ))
3381  zexp4 = exp_tbl(itind)
3382  endif
3383 
3384  ztdbt0 = zexp4 * ztdbt0
3385 
3386  else ! if_cldfmc_block --- it is a clear layer
3387 
3388 ! --- ... direct beam transmittance
3389  ztdbt(k) = zldbt(kp) * ztdbt(kp)
3390 
3391 ! --- ... pre-delta-scaling clear and cloudy direct beam transmittance
3392  ztdbt0 = zldbt0(k) * ztdbt0
3393 
3394  endif ! end if_cldfmc_block
3395  enddo ! end do_k_loop
3396 
3397 ! --- ... perform vertical quadrature
3398 
3399  call swflux &
3400 ! --- inputs:
3401  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3402  & nlay, nlp1, &
3403 ! --- outputs:
3404  & zfu, zfd &
3405  & )
3406 
3407 ! --- ... compute upward and downward fluxes at levels
3408  do k = 1, nlp1
3409  fxupc(k,ib) = fxupc(k,ib) + zsolar*zfu(k)
3410  fxdnc(k,ib) = fxdnc(k,ib) + zsolar*zfd(k)
3411  enddo
3412 
3413 !! --- ... surface downward beam/diffused flux components
3414  zb1 = zsolar*ztdbt0
3415  zb2 = zsolar*(zfd(1) - ztdbt0)
3416 
3417  if (ibd /= 0) then
3418  sfbmc(ibd) = sfbmc(ibd) + zb1
3419  sfdfc(ibd) = sfdfc(ibd) + zb2
3420  else
3421  zf1 = 0.5 * zb1
3422  zf2 = 0.5 * zb2
3423  sfbmc(1) = sfbmc(1) + zf1
3424  sfdfc(1) = sfdfc(1) + zf2
3425  sfbmc(2) = sfbmc(2) + zf1
3426  sfdfc(2) = sfdfc(2) + zf2
3427  endif
3428 ! sfbmc(ibd) = sfbmc(ibd) + zsolar*ztdbt0
3429 ! sfdfc(ibd) = sfdfc(ibd) + zsolar*(zfd(1) - ztdbt0)
3430 
3431  endif ! end if_cf1_block
3432 
3433  enddo lab_do_jg
3434 
3435 ! --- ... end of g-point loop
3436 
3437  do ib = 1, nbdsw
3438  ftoadc = ftoadc + fxdn0(nlp1,ib)
3439  ftoau0 = ftoau0 + fxup0(nlp1,ib)
3440  fsfcu0 = fsfcu0 + fxup0(1,ib)
3441  fsfcd0 = fsfcd0 + fxdn0(1,ib)
3442  enddo
3443 
3444 !! --- ... uv-b surface downward flux
3445  ibd = nuvb - nblow + 1
3446  suvbf0 = fxdn0(1,ibd)
3447 
3448  if ( cf1 <= eps ) then ! clear column, set total-sky=clear-sky fluxes
3449  do ib = 1, nbdsw
3450  do k = 1, nlp1
3451  fxupc(k,ib) = fxup0(k,ib)
3452  fxdnc(k,ib) = fxdn0(k,ib)
3453  enddo
3454  enddo
3455 
3456  ftoauc = ftoau0
3457  fsfcuc = fsfcu0
3458  fsfcdc = fsfcd0
3459 
3460 !! --- ... surface downward beam/diffused flux components
3461  sfbmc(1) = sfbm0(1)
3462  sfdfc(1) = sfdf0(1)
3463  sfbmc(2) = sfbm0(2)
3464  sfdfc(2) = sfdf0(2)
3465 
3466 !! --- ... uv-b surface downward flux
3467  suvbfc = suvbf0
3468  else ! cloudy column, compute total-sky fluxes
3469  do ib = 1, nbdsw
3470  ftoauc = ftoauc + fxupc(nlp1,ib)
3471  fsfcuc = fsfcuc + fxupc(1,ib)
3472  fsfcdc = fsfcdc + fxdnc(1,ib)
3473  enddo
3474 
3475 !! --- ... uv-b surface downward flux
3476  suvbfc = fxdnc(1,ibd)
3477  endif ! end if_cf1_block
3478 
3479  return
3480 !...................................

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine module_radsw_main::swflux ( )
private

Definition at line 3488 of file radsw_main.f.

References f_one, and f_zero.

Referenced by spcvrtc(), and spcvrtm().

3488 ! --- inputs:
3489  & ( zrefb,zrefd,ztrab,ztrad,zldbt,ztdbt, &
3490  & nlay, nlp1, &
3491 ! --- outputs:
3492  & zfu, zfd &
3493  & )
3494 
3495 ! =================== program usage description =================== !
3496 ! !
3497 ! purpose: computes the upward and downward radiation fluxes !
3498 ! !
3499 ! interface: "swflux" is called by "spcvrc" and "spcvrm" !
3500 ! !
3501 ! subroutines called : none !
3502 ! !
3503 ! ==================== defination of variables ==================== !
3504 ! !
3505 ! input variables: !
3506 ! zrefb(NLP1) - layer direct beam reflectivity !
3507 ! zrefd(NLP1) - layer diffuse reflectivity !
3508 ! ztrab(NLP1) - layer direct beam transmissivity !
3509 ! ztrad(NLP1) - layer diffuse transmissivity !
3510 ! zldbt(NLP1) - layer mean beam transmittance !
3511 ! ztdbt(NLP1) - total beam transmittance at levels !
3512 ! NLAY, NLP1 - number of layers/levels !
3513 ! !
3514 ! output variables: !
3515 ! zfu (NLP1) - upward flux at layer interface !
3516 ! zfd (NLP1) - downward flux at layer interface !
3517 ! !
3518 ! ******************************************************************* !
3519 ! ====================== end of description block ================= !
3520 
3521 ! --- inputs:
3522  integer, intent(in) :: nlay, nlp1
3523 
3524  real (kind=kind_phys), dimension(nlp1), intent(in) :: zrefb, &
3525  & zrefd, ztrab, ztrad, ztdbt, zldbt
3526 
3527 ! --- outputs:
3528  real (kind=kind_phys), dimension(nlp1), intent(out) :: zfu, zfd
3529 
3530 ! --- locals:
3531  real (kind=kind_phys), dimension(nlp1) :: zrupb,zrupd,zrdnd,ztdn
3532 
3533  real (kind=kind_phys) :: zden1
3534 
3535  integer :: k, kp
3536 !
3537 !===> ... begin here
3538 !
3539 
3540 ! --- ... link lowest layer with surface
3541 
3542  zrupb(1) = zrefb(1) ! direct beam
3543  zrupd(1) = zrefd(1) ! diffused
3544 
3545 ! --- ... pass from bottom to top
3546 
3547  do k = 1, nlay
3548  kp = k + 1
3549 
3550  zden1 = f_one / ( f_one - zrupd(k)*zrefd(kp) )
3551  zrupb(kp) = zrefb(kp) + ( ztrad(kp) * &
3552  & ( (ztrab(kp) - zldbt(kp))*zrupd(k) + &
3553  & zldbt(kp)*zrupb(k)) ) * zden1
3554  zrupd(kp) = zrefd(kp) + ztrad(kp)*ztrad(kp)*zrupd(k)*zden1
3555  enddo
3556 
3557 ! --- ... upper boundary conditions
3558 
3559  ztdn(nlp1) = f_one
3560  zrdnd(nlp1) = f_zero
3561  ztdn(nlay) = ztrab(nlp1)
3562  zrdnd(nlay) = zrefd(nlp1)
3563 
3564 ! --- ... pass from top to bottom
3565 
3566  do k = nlay, 2, -1
3567  zden1 = f_one / (f_one - zrefd(k)*zrdnd(k))
3568  ztdn(k-1) = ztdbt(k)*ztrab(k) + ( ztrad(k) * &
3569  & ( (ztdn(k) - ztdbt(k)) + ztdbt(k) * &
3570  & zrefb(k)*zrdnd(k) )) * zden1
3571  zrdnd(k-1) = zrefd(k) + ztrad(k)*ztrad(k)*zrdnd(k)*zden1
3572  enddo
3573 
3574 ! --- ... up and down-welling fluxes at levels
3575 
3576  do k = 1, nlp1
3577  zden1 = f_one / (f_one - zrdnd(k)*zrupd(k))
3578  zfu(k) = ( ztdbt(k)*zrupb(k) + &
3579  & (ztdn(k) - ztdbt(k))*zrupd(k) ) * zden1
3580  zfd(k) = ztdbt(k) + ( ztdn(k) - ztdbt(k) + &
3581  & ztdbt(k)*zrupb(k)*zrdnd(k) ) * zden1
3582  enddo
3583 
3584  return
3585 !...................................

Here is the caller graph for this function:

subroutine, public module_radsw_main::swrad ( )
Parameters
[in]plyr(npts,nlay),model layer mean pressure in mb
[in]plvl(npts,nlp1),model level pressure in mb
[in]tlyr(npts,nlay),model layer mean temperature in K
[in]tlvl(npts,nlp1),model level temperature in K (not in use)
[in]qlyr(npts,nlay),layer specific humidity in gm/gm *see inside
[in]olyr(npts,nlay),layer ozone concentration in gm/gm
[in]gasvmr(npts,nlay,:),atmospheric constent gases:(check module_radiation_gases for definition)
(:,:,1) - co2 volume mixing ratio
(:,:,2) - n2o volume mixing ratio
(:,:,3) - ch4 volume mixing ratio
(:,:,4) - o2 volume mixing ratio
(:,:,5) - co volume mixing ratio (not used)
(:,:,6) - cfc11 volume mixing ratio (not used)
(:,:,7) - cfc12 volume mixing ratio (not used)
(:,:,8) - cfc22 volume mixing ratio (not used)
(:,:,9) - ccl4 volume mixing ratio (not used)
[in]clouds(npts,nlay,:),cloud profile(check module_radiation_clouds for definition)
— for iswcliq > 0 —
(:,:,1) - layer total cloud fraction
(:,:,2) - layer in-cloud liq water path (g/m**2)
(:,:,3) - mean eff radius for liq cloud (micron)
(:,:,4) - layer in-cloud ice water path (g/m**2)
(:,:,5) - mean eff radius for ice cloud (micron)
(:,:,6) - layer rain drop water path (g/m**2)
(:,:,7) - mean eff radius for rain drop (micron)
(:,:,8) - layer snow flake water path (g/m**2)
(:,:,9) - mean eff radius for snow flake (micron)
— for iswcliq = 0 —
(:,:,1) - layer total cloud fraction
(:,:,2) - layer cloud optical depth
(:,:,3) - layer cloud single scattering albedo
(:,:,4) - layer cloud asymmetry factor
[in]icseed(npts),auxiliary special cloud related array.when module variable isubcsw=2, it provides permutation seed for each column profile that are used for generating random numbers.when isubcsw /=2, it will not be used.
[in]aerosols(npts,nlay,nbdsw,:),aerosol optical properties (check module_radiation_aerosols for definition)
(:,:,:,1) - optical depth
(:,:,:,2) - single scattering albedo
(:,:,:,3) - asymmetry parameter
[in]sfcalb(npts, : ),surface albedo in fraction(check module_radiation_surface for definition)
(:,1) - near ir direct beam albedo
(:,2) - near ir diffused albedo
(:,3) - uv+vis direct beam albedo
(:,4) - uv+vis diffused albedo
[in]cosz(npts),cosine of solar zenith angle
[in]solconsolar constant(w/m**2)
[in]NDAYnum of daytime points
[in]idxday(npts),index array for daytime points
[in]nptsnumber of horizontal points
[in]nlay,nlp1vertical layer/lavel numbers
[in]lprntlogical check print flag
[out]hswc(npts,nlay),total sky heating rates (k/sec or k/day)
[out]topflx(npts),radiation fluxes at toa (w/m**2), components (check module_radsw_parameters for definition)
upfxc - total sky upward flux at toa
dnflx - total sky downward flux at toa
upfx0 - clear sky upward flux at toa
[out]sfcflx(npts),radiation fluxes at sfc (w/m**2), components (check module_radsw_parameters for definition)
upfxc - total sky upward flux at sfc
dnfxc - total sky downward flux at sfc
upfx0 - clear sky upward flux at sfc
dnfx0 - clear sky downward flux at sfc
Optional:
[out]hswb(npts,nlay,nbdsw),spectral band total sky heating rates
[out]hsw0(npts,nlay),clear sky heating rates (k/sec or k/day)
[out]flxprf(npts,nlp1),level radiation fluxes (w/m**2), components (check module_radsw_parameters for definition)
dnfxc - total sky downward flux at interface
upfxc - total sky upward flux at interface
dnfx0 - clear sky downward flux at interface
upfx0 - clear sky upward flux at interface
[out]fdncmp(npts),component surface downward fluxes (w/m**2)(check module_radsw_parameters for definition) !
uvbfc - total sky downward uv-b flux at sfc
uvbf0 - clear sky downward uv-b flux at sfc
nirbm - downward surface nir direct beam flux
nirdf - downward surface nir diffused flux
visbm - downward surface uv+vis direct beam flux
visdf - downward surface uv+vis diffused flux

General Algorithm

  1. Compute solar constant adjustment factor (s0fac) according to solcon.
    \( s0fac = solcon/s0 \)
    where s0, the solar constant at toa in \( w/m^2 \), is hard-coded with each spectran band, the total flux is about 1368.22 \( w/m^2 \).
  2. Initial output arrays (and optional) as zero
  3. Change random number seed value for each radiation invocation (isubcsw =1 or 2).
  4. Prepare surface albedo: bm,df - dir,dif; 1,2 - nir,uvv
  5. Prepare atmospheric profile for use in rrtm:
    • pavel(nlay):layer pressure (mb)
    • tavel(nlay):layer temperature (k)
    • delp (nlay):layer pressure thickness (mb)
  6. Set absorber and gas column amount, convert from volume mixing ratio to molec/cm2 based on coldry (scaled to 1.0e-20)
    • colamt(nlay,maxgas):column amounts of absorbing gases 1 to maxgas are for h2o,co2,o3,n2o,ch4,o2,co, respectively ( \( mol/cm^2 \))
  7. Read aerosol optical properties from aerosols (:,:,:,:)
  8. Read cloud optical properties from clouds(:,:,:)
  9. Compute fractions of clear sky view in different overlapping scenario:
    • random overlapping
    • max/ran overlapping
    • maximum overlapping
  10. For cloudy sky column, call cldprop (Cloud Optical Property in SW) to Compute the cloud optical properties for each cloudy layer
  11. Calling setcoef (SW Coefficients calculation)to compute various coefficients needed in radiative transfer calculations.
  12. Calling taumol to calculate optical depths for gaseous absorption and rayleigh scattering
    - SUBROUTINE taumol calls the subroutines taugbn (where n goes from 16 to 29). taugbn calculates the optical depths and Planck fractions per g-value and layer for band n.
  13. call the 2-stream radiation transfer model
  14. Save outputs

Definition at line 440 of file radsw_main.f.

References amdo3, amdw, cldprop(), physcons::con_amd, physcons::con_amw, physcons::con_avgd, physcons::con_g, f_one, f_zero, ftiny, heatfac, physparam::iovrsw, ipsdsw0, physparam::isubcsw, physparam::iswcliq, physparam::iswrgas, physparam::ivflip, lfdncmp, lflxprf, lhsw0, lhswb, module_radsw_parameters::nbdsw, oneminus, s0, setcoef(), spcvrtc(), spcvrtm(), and taumol().

Referenced by module_radiation_driver::grrad().

440 
441 ! --- inputs:
442  & ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
443  & clouds,icseed,aerosols,sfcalb, &
444  & cosz,solcon,nday,idxday, &
445  & npts, nlay, nlp1, lprnt, &
446 ! --- outputs:
447  & hswc,topflx,sfcflx &
448 !! --- optional:
449  &, hsw0,hswb,flxprf,fdncmp &
450  & )
451 
452 ! ==================== defination of variables ==================== !
453 ! !
454 ! input variables: !
455 ! plyr (npts,nlay) : model layer mean pressure in mb !
456 ! plvl (npts,nlp1) : model level pressure in mb !
457 ! tlyr (npts,nlay) : model layer mean temperature in k !
458 ! tlvl (npts,nlp1) : model level temperature in k (not in use) !
459 ! qlyr (npts,nlay) : layer specific humidity in gm/gm *see inside !
460 ! olyr (npts,nlay) : layer ozone concentration in gm/gm !
461 ! gasvmr(npts,nlay,:): atmospheric constent gases: !
462 ! (check module_radiation_gases for definition) !
463 ! gasvmr(:,:,1) - co2 volume mixing ratio !
464 ! gasvmr(:,:,2) - n2o volume mixing ratio !
465 ! gasvmr(:,:,3) - ch4 volume mixing ratio !
466 ! gasvmr(:,:,4) - o2 volume mixing ratio !
467 ! gasvmr(:,:,5) - co volume mixing ratio (not used) !
468 ! gasvmr(:,:,6) - cfc11 volume mixing ratio (not used) !
469 ! gasvmr(:,:,7) - cfc12 volume mixing ratio (not used) !
470 ! gasvmr(:,:,8) - cfc22 volume mixing ratio (not used) !
471 ! gasvmr(:,:,9) - ccl4 volume mixing ratio (not used) !
472 ! clouds(npts,nlay,:): cloud profile !
473 ! (check module_radiation_clouds for definition) !
474 ! --- for iswcliq > 0 --- !
475 ! clouds(:,:,1) - layer total cloud fraction !
476 ! clouds(:,:,2) - layer in-cloud liq water path (g/m**2) !
477 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
478 ! clouds(:,:,4) - layer in-cloud ice water path (g/m**2) !
479 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
480 ! clouds(:,:,6) - layer rain drop water path (g/m**2) !
481 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
482 ! clouds(:,:,8) - layer snow flake water path (g/m**2) !
483 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
484 ! --- for iswcliq = 0 --- !
485 ! clouds(:,:,1) - layer total cloud fraction !
486 ! clouds(:,:,2) - layer cloud optical depth !
487 ! clouds(:,:,3) - layer cloud single scattering albedo !
488 ! clouds(:,:,4) - layer cloud asymmetry factor !
489 ! icseed(npts) : auxiliary special cloud related array !
490 ! when module variable isubcsw=2, it provides !
491 ! permutation seed for each column profile that !
492 ! are used for generating random numbers. !
493 ! when isubcsw /=2, it will not be used. !
494 ! aerosols(npts,nlay,nbdsw,:) : aerosol optical properties !
495 ! (check module_radiation_aerosols for definition) !
496 ! (:,:,:,1) - optical depth !
497 ! (:,:,:,2) - single scattering albedo !
498 ! (:,:,:,3) - asymmetry parameter !
499 ! sfcalb(npts, : ) : surface albedo in fraction !
500 ! (check module_radiation_surface for definition) !
501 ! ( :, 1 ) - near ir direct beam albedo !
502 ! ( :, 2 ) - near ir diffused albedo !
503 ! ( :, 3 ) - uv+vis direct beam albedo !
504 ! ( :, 4 ) - uv+vis diffused albedo !
505 ! cosz (npts) : cosine of solar zenith angle !
506 ! solcon : solar constant (w/m**2) !
507 ! NDAY : num of daytime points !
508 ! idxday(npts) : index array for daytime points !
509 ! npts : number of horizontal points !
510 ! nlay,nlp1 : vertical layer/lavel numbers !
511 ! lprnt : logical check print flag !
512 ! !
513 ! output variables: !
514 ! hswc (npts,nlay): total sky heating rates (k/sec or k/day) !
515 ! topflx(npts) : radiation fluxes at toa (w/m**2), components: !
516 ! (check module_radsw_parameters for definition) !
517 ! upfxc - total sky upward flux at toa !
518 ! dnflx - total sky downward flux at toa !
519 ! upfx0 - clear sky upward flux at toa !
520 ! sfcflx(npts) : radiation fluxes at sfc (w/m**2), components: !
521 ! (check module_radsw_parameters for definition) !
522 ! upfxc - total sky upward flux at sfc !
523 ! dnfxc - total sky downward flux at sfc !
524 ! upfx0 - clear sky upward flux at sfc !
525 ! dnfx0 - clear sky downward flux at sfc !
526 ! !
527 !!optional outputs variables: !
528 ! hswb(npts,nlay,nbdsw): spectral band total sky heating rates !
529 ! hsw0 (npts,nlay): clear sky heating rates (k/sec or k/day) !
530 ! flxprf(npts,nlp1): level radiation fluxes (w/m**2), components: !
531 ! (check module_radsw_parameters for definition) !
532 ! dnfxc - total sky downward flux at interface !
533 ! upfxc - total sky upward flux at interface !
534 ! dnfx0 - clear sky downward flux at interface !
535 ! upfx0 - clear sky upward flux at interface !
536 ! fdncmp(npts) : component surface downward fluxes (w/m**2): !
537 ! (check module_radsw_parameters for definition) !
538 ! uvbfc - total sky downward uv-b flux at sfc !
539 ! uvbf0 - clear sky downward uv-b flux at sfc !
540 ! nirbm - downward surface nir direct beam flux !
541 ! nirdf - downward surface nir diffused flux !
542 ! visbm - downward surface uv+vis direct beam flux !
543 ! visdf - downward surface uv+vis diffused flux !
544 ! !
545 ! external module variables: (in physparam) !
546 ! iswrgas - control flag for rare gases (ch4,n2o,o2, etc.) !
547 ! =0: do not include rare gases !
548 ! >0: include all rare gases !
549 ! iswcliq - control flag for liq-cloud optical properties !
550 ! =0: input cloud optical depth, fixed ssa, asy !
551 ! =1: use hu and stamnes(1993) method for liq cld !
552 ! =2: not used !
553 ! iswcice - control flag for ice-cloud optical properties !
554 ! *** if iswcliq==0, iswcice is ignored !
555 ! =1: use ebert and curry (1992) scheme for ice clouds !
556 ! =2: use streamer v3.0 (2001) method for ice clouds !
557 ! =3: use fu's method (1996) for ice clouds !
558 ! iswmode - control flag for 2-stream transfer scheme !
559 ! =1; delta-eddington (joseph et al., 1976) !
560 ! =2: pifm (zdunkowski et al., 1980) !
561 ! =3: discrete ordinates (liou, 1973) !
562 ! isubcsw - sub-column cloud approximation control flag !
563 ! =0: no sub-col cld treatment, use grid-mean cld quantities !
564 ! =1: mcica sub-col, prescribed seeds to get random numbers !
565 ! =2: mcica sub-col, providing array icseed for random numbers!
566 ! iovrsw - cloud overlapping control flag !
567 ! =0: random overlapping clouds !
568 ! =1: maximum/random overlapping clouds !
569 ! =2: maximum overlap cloud !
570 ! ivflip - control flg for direction of vertical index !
571 ! =0: index from toa to surface !
572 ! =1: index from surface to toa !
573 ! !
574 ! module parameters, control variables: !
575 ! nblow,nbhgh - lower and upper limits of spectral bands !
576 ! maxgas - maximum number of absorbing gaseous !
577 ! ngptsw - total number of g-point subintervals !
578 ! ng## - number of g-points in band (##=16-29) !
579 ! ngb(ngptsw) - band indices for each g-point !
580 ! bpade - pade approximation constant (1/0.278) !
581 ! nspa,nspb(nblow:nbhgh) !
582 ! - number of lower/upper ref atm's per band !
583 ! ipsdsw0 - permutation seed for mcica sub-col clds !
584 ! !
585 ! major local variables: !
586 ! pavel (nlay) - layer pressures (mb) !
587 ! delp (nlay) - layer pressure thickness (mb) !
588 ! tavel (nlay) - layer temperatures (k) !
589 ! coldry (nlay) - dry air column amount !
590 ! (1.e-20*molecules/cm**2) !
591 ! cldfrc (nlay) - layer cloud fraction (norm by tot cld) !
592 ! cldfmc (nlay,ngptsw) - layer cloud fraction for g-point !
593 ! taucw (nlay,nbdsw) - cloud optical depth !
594 ! ssacw (nlay,nbdsw) - cloud single scattering albedo (weighted) !
595 ! asycw (nlay,nbdsw) - cloud asymmetry factor (weighted) !
596 ! tauaer (nlay,nbdsw) - aerosol optical depths !
597 ! ssaaer (nlay,nbdsw) - aerosol single scattering albedo !
598 ! asyaer (nlay,nbdsw) - aerosol asymmetry factor !
599 ! colamt (nlay,maxgas) - column amounts of absorbing gases !
600 ! 1 to maxgas are for h2o, co2, o3, n2o, !
601 ! ch4, o2, co, respectively (mol/cm**2) !
602 ! facij (nlay) - indicator of interpolation factors !
603 ! =0/1: indicate lower/higher temp & height !
604 ! selffac(nlay) - scale factor for self-continuum, equals !
605 ! (w.v. density)/(atm density at 296K,1013 mb) !
606 ! selffrac(nlay) - factor for temp interpolation of ref !
607 ! self-continuum data !
608 ! indself(nlay) - index of the lower two appropriate ref !
609 ! temp for the self-continuum interpolation !
610 ! forfac (nlay) - scale factor for w.v. foreign-continuum !
611 ! forfrac(nlay) - factor for temp interpolation of ref !
612 ! w.v. foreign-continuum data !
613 ! indfor (nlay) - index of the lower two appropriate ref !
614 ! temp for the foreign-continuum interp !
615 ! laytrop - layer at which switch is made from one !
616 ! combination of key species to another !
617 ! jp(nlay),jt(nlay),jt1(nlay) !
618 ! - lookup table indexes !
619 ! flxucb(nlp1,nbdsw) - spectral bnd total-sky upward flx (w/m2) !
620 ! flxdcb(nlp1,nbdsw) - spectral bnd total-sky downward flx (w/m2)!
621 ! flxu0b(nlp1,nbdsw) - spectral bnd clear-sky upward flx (w/m2) !
622 ! flxd0b(nlp1,nbdsw) - spectral b d clear-sky downward flx (w/m2)!
623 ! !
624 ! !
625 ! ===================== end of definitions ==================== !
626 
627 ! --- inputs:
628  integer, intent(in) :: npts, nlay, nlp1, nday
629 
630  integer, dimension(:), intent(in) :: idxday, icseed
631 
632  logical, intent(in) :: lprnt
633 
634  real (kind=kind_phys), dimension(npts,nlp1), intent(in) :: &
635  & plvl, tlvl
636  real (kind=kind_phys), dimension(npts,nlay), intent(in) :: &
637  & plyr, tlyr, qlyr, olyr
638  real (kind=kind_phys), dimension(npts,4), intent(in) :: sfcalb
639 
640  real (kind=kind_phys), dimension(npts,nlay,9),intent(in):: gasvmr
641  real (kind=kind_phys), dimension(npts,nlay,9),intent(in):: clouds
642  real (kind=kind_phys), dimension(npts,nlay,nbdsw,3),intent(in):: &
643  & aerosols
644 
645  real (kind=kind_phys), intent(in) :: cosz(npts), solcon
646 
647 ! --- outputs:
648  real (kind=kind_phys), dimension(npts,nlay), intent(out) :: hswc
649 
650  type(topfsw_type), dimension(npts), intent(out) :: topflx
651  type(sfcfsw_type), dimension(npts), intent(out) :: sfcflx
652 
653 !! --- optional outputs:
654  real (kind=kind_phys), dimension(npts,nlay,nbdsw), optional, &
655  & intent(out) :: hswb
656 
657  real (kind=kind_phys), dimension(npts,nlay), optional, &
658  & intent(out) :: hsw0
659  type (profsw_type), dimension(npts,nlp1), optional, &
660  & intent(out) :: flxprf
661  type (cmpfsw_type), dimension(npts), optional, &
662  & intent(out) :: fdncmp
663 
664 ! --- locals:
665  real (kind=kind_phys), dimension(nlay,ngptsw) :: cldfmc, &
666  & taug, taur
667  real (kind=kind_phys), dimension(nlp1,nbdsw):: fxupc, fxdnc, &
668  & fxup0, fxdn0
669 
670  real (kind=kind_phys), dimension(nlay,nbdsw) :: &
671  & tauae, ssaae, asyae, taucw, ssacw, asycw
672 
673  real (kind=kind_phys), dimension(ngptsw) :: sfluxzen
674 
675  real (kind=kind_phys), dimension(nlay) :: cldfrc, delp, &
676  & pavel, tavel, coldry, colmol, h2ovmr, o3vmr, temcol, &
677  & cliqp, reliq, cicep, reice, cdat1, cdat2, cdat3, cdat4, &
678  & cfrac, fac00, fac01, fac10, fac11, forfac, forfrac, &
679  & selffac, selffrac, rfdelp
680 
681  real (kind=kind_phys), dimension(nlp1) :: fnet, flxdc, flxuc, &
682  & flxd0, flxu0
683 
684  real (kind=kind_phys), dimension(2) :: albbm, albdf, sfbmc, &
685  & sfbm0, sfdfc, sfdf0
686 
687  real (kind=kind_phys) :: cosz1, sntz1, tem0, tem1, tem2, s0fac, &
688  & ssolar, zcf0, zcf1, ftoau0, ftoauc, ftoadc, &
689  & fsfcu0, fsfcuc, fsfcd0, fsfcdc, suvbfc, suvbf0
690 
691 ! --- column amount of absorbing gases:
692 ! (:,m) m = 1-h2o, 2-co2, 3-o3, 4-n2o, 5-ch4, 6-o2, 7-co
693  real (kind=kind_phys) :: colamt(nlay,maxgas)
694 
695  integer, dimension(npts) :: ipseed
696  integer, dimension(nlay) :: indfor, indself, jp, jt, jt1
697 
698  integer :: i, ib, ipt, j1, k, kk, laytrop, mb
699 !
700 !===> ... begin here
701 !
702 
703  lhswb = present ( hswb )
704  lhsw0 = present ( hsw0 )
705  lflxprf= present ( flxprf )
706  lfdncmp= present ( fdncmp )
707 
712 ! --- ... compute solar constant adjustment factor according to solcon.
713 ! *** s0, the solar constant at toa in w/m**2, is hard-coded with
714 ! each spectra band, the total flux is about 1368.22 w/m**2.
715 
716  s0fac = solcon / s0
717 
719 ! --- ... initial output arrays
720 
721  hswc(:,:) = f_zero
722  topflx = topfsw_type( f_zero, f_zero, f_zero )
723  sfcflx = sfcfsw_type( f_zero, f_zero, f_zero, f_zero )
724 
725 !! --- ... initial optional outputs
726  if ( lflxprf ) then
727  flxprf = profsw_type( f_zero, f_zero, f_zero, f_zero )
728  endif
729 
730  if ( lfdncmp ) then
731  fdncmp = cmpfsw_type(f_zero,f_zero,f_zero,f_zero,f_zero,f_zero)
732  endif
733 
734  if ( lhsw0 ) then
735  hsw0(:,:) = f_zero
736  endif
737 
738  if ( lhswb ) then
739  hswb(:,:,:) = f_zero
740  endif
741 
743 ! --- ... change random number seed value for each radiation invocation
744 
745  if ( isubcsw == 1 ) then ! advance prescribed permutation seed
746  do i = 1, npts
747  ipseed(i) = ipsdsw0 + i
748  enddo
749  elseif ( isubcsw == 2 ) then ! use input array of permutaion seeds
750  do i = 1, npts
751  ipseed(i) = icseed(i)
752  enddo
753  endif
754 
755  if ( lprnt ) then
756  write(0,*)' In radsw, isubcsw, ipsdsw0,ipseed =', &
757  & isubcsw, ipsdsw0, ipseed
758  endif
759 
760 ! --- ... loop over each daytime grid point
761 
762  lab_do_ipt : do ipt = 1, nday
763 
764  j1 = idxday(ipt)
765 
766  cosz1 = cosz(j1)
767  sntz1 = f_one / cosz(j1)
768  ssolar = s0fac * cosz(j1)
769 
771 ! --- ... surface albedo: bm,df - dir,dif; 1,2 - nir,uvv
772  albbm(1) = sfcalb(j1,1)
773  albdf(1) = sfcalb(j1,2)
774  albbm(2) = sfcalb(j1,3)
775  albdf(2) = sfcalb(j1,4)
776 
783 
784 ! --- ... prepare atmospheric profile for use in rrtm
785 ! the vertical index of internal array is from surface to top
786 
787  if (ivflip == 0) then ! input from toa to sfc
788 
789  tem1 = 100.0 * con_g
790  tem2 = 1.0e-20 * 1.0e3 * con_avgd
791 
792  do k = 1, nlay
793  kk = nlp1 - k
794  pavel(k) = plyr(j1,kk)
795  tavel(k) = tlyr(j1,kk)
796  delp(k) = plvl(j1,kk+1) - plvl(j1,kk)
797 
798 ! --- ... set absorber amount
799 !test use
800 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw) ! input mass mixing ratio
801 ! h2ovmr(k)= max(f_zero,qlyr(j1,kk)) ! input vol mixing ratio
802 ! o3vmr (k)= max(f_zero,olyr(j1,kk)) ! input vol mixing ratio
803 !ncep model use
804  h2ovmr(k)= max(f_zero,qlyr(j1,kk)*amdw/(f_one-qlyr(j1,kk))) ! input specific humidity
805  o3vmr(k)= max(f_zero,olyr(j1,kk)*amdo3) ! input mass mixing ratio
806 
807  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
808  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
809  temcol(k) = 1.0e-12 * coldry(k)
810 
811  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
812  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,kk,1)) ! co2
813  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
814  colmol(k) = coldry(k) + colamt(k,1)
815  enddo
816 
817 ! --- ... set up gas column amount, convert from volume mixing ratio
818 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
819 
820  if (iswrgas > 0) then
821  do k = 1, nlay
822  kk = nlp1 - k
823  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,kk,2)) ! n2o
824  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,kk,3)) ! ch4
825  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,kk,4)) ! o2
826 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,kk,5)) ! co - notused
827  enddo
828  else
829  do k = 1, nlay
830  colamt(k,4) = temcol(k) ! n2o
831  colamt(k,5) = temcol(k) ! ch4
832  colamt(k,6) = temcol(k) ! o2
833 ! colamt(k,7) = temcol(k) ! co - notused
834  enddo
835  endif
836 
838 ! --- ... set aerosol optical properties
839 
840  do k = 1, nlay
841  kk = nlp1 - k
842  do ib = 1, nbdsw
843  tauae(k,ib) = aerosols(j1,kk,ib,1)
844  ssaae(k,ib) = aerosols(j1,kk,ib,2)
845  asyae(k,ib) = aerosols(j1,kk,ib,3)
846  enddo
847  enddo
848 
850  if (iswcliq > 0) then ! use prognostic cloud method
851  do k = 1, nlay
852  kk = nlp1 - k
853  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
854  cliqp(k) = clouds(j1,kk,2) ! cloud liq path
855  reliq(k) = clouds(j1,kk,3) ! liq partical effctive radius
856  cicep(k) = clouds(j1,kk,4) ! cloud ice path
857  reice(k) = clouds(j1,kk,5) ! ice partical effctive radius
858  cdat1(k) = clouds(j1,kk,6) ! cloud rain drop path
859  cdat2(k) = clouds(j1,kk,7) ! rain partical effctive radius
860  cdat3(k) = clouds(j1,kk,8) ! cloud snow path
861  cdat4(k) = clouds(j1,kk,9) ! snow partical effctive radius
862  enddo
863  else ! use diagnostic cloud method
864  do k = 1, nlay
865  kk = nlp1 - k
866  cfrac(k) = clouds(j1,kk,1) ! cloud fraction
867  cdat1(k) = clouds(j1,kk,2) ! cloud optical depth
868  cdat2(k) = clouds(j1,kk,3) ! cloud single scattering albedo
869  cdat3(k) = clouds(j1,kk,4) ! cloud asymmetry factor
870  enddo
871  endif ! end if_iswcliq
872 
873  else ! input from sfc to toa
874 
875  tem1 = 100.0 * con_g
876  tem2 = 1.0e-20 * 1.0e3 * con_avgd
877 
878  do k = 1, nlay
879  pavel(k) = plyr(j1,k)
880  tavel(k) = tlyr(j1,k)
881  delp(k) = plvl(j1,k) - plvl(j1,k+1)
882 
883 ! --- ... set absorber amount
884 !test use
885 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw) ! input mass mixing ratio
886 ! h2ovmr(k)= max(f_zero,qlyr(j1,k)) ! input vol mixing ratio
887 ! o3vmr (k)= max(f_zero,olyr(j1,k)) ! input vol mixing ratio
888 !ncep model use
889  h2ovmr(k)= max(f_zero,qlyr(j1,k)*amdw/(f_one-qlyr(j1,k))) ! input specific humidity
890  o3vmr(k)= max(f_zero,olyr(j1,k)*amdo3) ! input mass mixing ratio
891 
892  tem0 = (f_one - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
893  coldry(k) = tem2 * delp(k) / (tem1*tem0*(f_one + h2ovmr(k)))
894  temcol(k) = 1.0e-12 * coldry(k)
895 
896  colamt(k,1) = max(f_zero, coldry(k)*h2ovmr(k)) ! h2o
897  colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(j1,k,1)) ! co2
898  colamt(k,3) = max(f_zero, coldry(k)*o3vmr(k)) ! o3
899  colmol(k) = coldry(k) + colamt(k,1)
900  enddo
901 
902 
903  if (lprnt) then
904  if (ipt == 1) then
905  write(0,*)' pavel=',pavel
906  write(0,*)' tavel=',tavel
907  write(0,*)' delp=',delp
908  write(0,*)' h2ovmr=',h2ovmr*1000
909  write(0,*)' o3vmr=',o3vmr*1000000
910  endif
911  endif
912 
913 ! --- ... set up gas column amount, convert from volume mixing ratio
914 ! to molec/cm2 based on coldry (scaled to 1.0e-20)
915 
916  if (iswrgas > 0) then
917  do k = 1, nlay
918  colamt(k,4) = max(temcol(k), coldry(k)*gasvmr(j1,k,2)) ! n2o
919  colamt(k,5) = max(temcol(k), coldry(k)*gasvmr(j1,k,3)) ! ch4
920  colamt(k,6) = max(temcol(k), coldry(k)*gasvmr(j1,k,4)) ! o2
921 ! colamt(k,7) = max(temcol(k), coldry(k)*gasvmr(j1,k,5)) ! co - notused
922  enddo
923  else
924  do k = 1, nlay
925  colamt(k,4) = temcol(k) ! n2o
926  colamt(k,5) = temcol(k) ! ch4
927  colamt(k,6) = temcol(k) ! o2
928 ! colamt(k,7) = temcol(k) ! co - notused
929  enddo
930  endif
931 
932 ! --- ... set aerosol optical properties
933 
934  do ib = 1, nbdsw
935  do k = 1, nlay
936  tauae(k,ib) = aerosols(j1,k,ib,1)
937  ssaae(k,ib) = aerosols(j1,k,ib,2)
938  asyae(k,ib) = aerosols(j1,k,ib,3)
939  enddo
940  enddo
941 
942  if (iswcliq > 0) then ! use prognostic cloud method
943  do k = 1, nlay
944  cfrac(k) = clouds(j1,k,1) ! cloud fraction
945  cliqp(k) = clouds(j1,k,2) ! cloud liq path
946  reliq(k) = clouds(j1,k,3) ! liq partical effctive radius
947  cicep(k) = clouds(j1,k,4) ! cloud ice path
948  reice(k) = clouds(j1,k,5) ! ice partical effctive radius
949  cdat1(k) = clouds(j1,k,6) ! cloud rain drop path
950  cdat2(k) = clouds(j1,k,7) ! rain partical effctive radius
951  cdat3(k) = clouds(j1,k,8) ! cloud snow path
952  cdat4(k) = clouds(j1,k,9) ! snow partical effctive radius
953  enddo
954  else ! use diagnostic cloud method
955  do k = 1, nlay
956  cfrac(k) = clouds(j1,k,1) ! cloud fraction
957  cdat1(k) = clouds(j1,k,2) ! cloud optical depth
958  cdat2(k) = clouds(j1,k,3) ! cloud single scattering albedo
959  cdat3(k) = clouds(j1,k,4) ! cloud asymmetry factor
960  enddo
961  endif ! end if_iswcliq
962 
963  endif ! if_ivflip
964 
969 ! --- ... compute fractions of clear sky view
970 
971  zcf0 = f_one
972  zcf1 = f_one
973  if (iovrsw == 0) then ! random overlapping
974  do k = 1, nlay
975  zcf0 = zcf0 * (f_one - cfrac(k))
976  enddo
977  else if (iovrsw == 1) then ! max/ran overlapping
978  do k = 1, nlay
979  if (cfrac(k) > ftiny) then ! cloudy layer
980  zcf1 = min( zcf1, f_one-cfrac(k) )
981  elseif (zcf1 < f_one) then ! clear layer
982  zcf0 = zcf0 * zcf1
983  zcf1 = f_one
984  endif
985  enddo
986  zcf0 = zcf0 * zcf1
987  else if (iovrsw == 2) then ! maximum overlapping
988  do k = 1, nlay
989  zcf0 = min( zcf0, f_one-cfrac(k) )
990  enddo
991  endif
992 
993  if (zcf0 <= ftiny) zcf0 = f_zero
994  if (zcf0 > oneminus) zcf0 = f_one
995  zcf1 = f_one - zcf0
996 
998 
999 ! --- ... compute cloud optical properties
1000 
1001  if (zcf1 > f_zero) then ! cloudy sky column
1002 
1003  call cldprop &
1004 ! --- inputs:
1005  & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, &
1006  & zcf1, nlay, ipseed(j1), &
1007 ! --- outputs:
1008  & taucw, ssacw, asycw, cldfrc, cldfmc &
1009  & )
1010 
1011  else ! clear sky column
1012  cldfrc(:) = f_zero
1013  cldfmc(:,:)= f_zero
1014  do i = 1, nbdsw
1015  do k = 1, nlay
1016  taucw(k,i) = f_zero
1017  ssacw(k,i) = f_zero
1018  asycw(k,i) = f_zero
1019  enddo
1020  enddo
1021  endif ! end if_zcf1_block
1022 
1024  call setcoef &
1025 ! --- inputs:
1026  & ( pavel,tavel,h2ovmr, nlay,nlp1, &
1027 ! --- outputs:
1028  & laytrop,jp,jt,jt1,fac00,fac01,fac10,fac11, &
1029  & selffac,selffrac,indself,forfac,forfrac,indfor &
1030  & )
1031 
1036 ! --- ... calculate optical depths for gaseous absorption and Rayleigh
1037 ! scattering
1038 
1039  call taumol &
1040 ! --- inputs:
1041  & ( colamt,colmol,fac00,fac01,fac10,fac11,jp,jt,jt1,laytrop, &
1042  & forfac,forfrac,indfor,selffac,selffrac,indself, nlay, &
1043 ! --- outputs:
1044  & sfluxzen, taug, taur &
1045  & )
1046 
1050 ! --- ... call the 2-stream radiation transfer model
1051 
1052  if ( isubcsw <= 0 ) then ! use standard cloud scheme
1053 
1054  call spcvrtc &
1055 ! --- inputs:
1056  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfrc, &
1057  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1058  & nlay, nlp1, &
1059 ! --- outputs:
1060  & fxupc,fxdnc,fxup0,fxdn0, &
1061  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1062  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1063  & )
1064 
1065  else ! use mcica cloud scheme
1066 
1067  call spcvrtm &
1068 ! --- inputs:
1069  & ( ssolar,cosz1,sntz1,albbm,albdf,sfluxzen,cldfmc, &
1070  & zcf1,zcf0,taug,taur,tauae,ssaae,asyae,taucw,ssacw,asycw, &
1071  & nlay, nlp1, &
1072 ! --- outputs:
1073  & fxupc,fxdnc,fxup0,fxdn0, &
1074  & ftoauc,ftoau0,ftoadc,fsfcuc,fsfcu0,fsfcdc,fsfcd0, &
1075  & sfbmc,sfdfc,sfbm0,sfdf0,suvbfc,suvbf0 &
1076  & )
1077 
1078  endif
1080 ! --- ... sum up total spectral fluxes for total-sky
1081 
1082  do k = 1, nlp1
1083  flxuc(k) = f_zero
1084  flxdc(k) = f_zero
1085 
1086  do ib = 1, nbdsw
1087  flxuc(k) = flxuc(k) + fxupc(k,ib)
1088  flxdc(k) = flxdc(k) + fxdnc(k,ib)
1089  enddo
1090  enddo
1091 
1092 !! --- ... optional clear sky fluxes
1093 
1094  if ( lhsw0 .or. lflxprf ) then
1095  do k = 1, nlp1
1096  flxu0(k) = f_zero
1097  flxd0(k) = f_zero
1098 
1099  do ib = 1, nbdsw
1100  flxu0(k) = flxu0(k) + fxup0(k,ib)
1101  flxd0(k) = flxd0(k) + fxdn0(k,ib)
1102  enddo
1103  enddo
1104  endif
1105 
1106 ! --- ... prepare for final outputs
1107 
1108  do k = 1, nlay
1109  rfdelp(k) = heatfac / delp(k)
1110  enddo
1111 
1112  if ( lfdncmp ) then
1113 !! --- ... optional uv-b surface downward flux
1114  fdncmp(j1)%uvbf0 = suvbf0
1115  fdncmp(j1)%uvbfc = suvbfc
1116 
1117 !! --- ... optional beam and diffuse sfc fluxes
1118  fdncmp(j1)%nirbm = sfbmc(1)
1119  fdncmp(j1)%nirdf = sfdfc(1)
1120  fdncmp(j1)%visbm = sfbmc(2)
1121  fdncmp(j1)%visdf = sfdfc(2)
1122  endif ! end if_lfdncmp
1123 
1124 ! --- ... toa and sfc fluxes
1125 
1126  topflx(j1)%upfxc = ftoauc
1127  topflx(j1)%dnfxc = ftoadc
1128  topflx(j1)%upfx0 = ftoau0
1129 
1130  sfcflx(j1)%upfxc = fsfcuc
1131  sfcflx(j1)%dnfxc = fsfcdc
1132  sfcflx(j1)%upfx0 = fsfcu0
1133  sfcflx(j1)%dnfx0 = fsfcd0
1134 
1135  if (ivflip == 0) then ! output from toa to sfc
1136 
1137 ! --- ... compute heating rates
1138 
1139  fnet(1) = flxdc(1) - flxuc(1)
1140 
1141  do k = 2, nlp1
1142  kk = nlp1 - k + 1
1143  fnet(k) = flxdc(k) - flxuc(k)
1144  hswc(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1145  enddo
1146 
1147 !! --- ... optional flux profiles
1148 
1149  if ( lflxprf ) then
1150  do k = 1, nlp1
1151  kk = nlp1 - k + 1
1152  flxprf(j1,kk)%upfxc = flxuc(k)
1153  flxprf(j1,kk)%dnfxc = flxdc(k)
1154  flxprf(j1,kk)%upfx0 = flxu0(k)
1155  flxprf(j1,kk)%dnfx0 = flxd0(k)
1156  enddo
1157  endif
1158 
1159 !! --- ... optional clear sky heating rates
1160 
1161  if ( lhsw0 ) then
1162  fnet(1) = flxd0(1) - flxu0(1)
1163 
1164  do k = 2, nlp1
1165  kk = nlp1 - k + 1
1166  fnet(k) = flxd0(k) - flxu0(k)
1167  hsw0(j1,kk) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1168  enddo
1169  endif
1170 
1171 !! --- ... optional spectral band heating rates
1172 
1173  if ( lhswb ) then
1174  do mb = 1, nbdsw
1175  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1176 
1177  do k = 2, nlp1
1178  kk = nlp1 - k + 1
1179  fnet(k) = fxdnc(k,mb) - fxupc(k,mb)
1180  hswb(j1,kk,mb) = (fnet(k) - fnet(k-1)) * rfdelp(k-1)
1181  enddo
1182  enddo
1183  endif
1184 
1185  else ! output from sfc to toa
1186 
1187 ! --- ... compute heating rates
1188 
1189  fnet(1) = flxdc(1) - flxuc(1)
1190 
1191  do k = 2, nlp1
1192  fnet(k) = flxdc(k) - flxuc(k)
1193  hswc(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1194  enddo
1195 
1196 !! --- ... optional flux profiles
1197 
1198  if ( lflxprf ) then
1199  do k = 1, nlp1
1200  flxprf(j1,k)%upfxc = flxuc(k)
1201  flxprf(j1,k)%dnfxc = flxdc(k)
1202  flxprf(j1,k)%upfx0 = flxu0(k)
1203  flxprf(j1,k)%dnfx0 = flxd0(k)
1204  enddo
1205  endif
1206 
1207 !! --- ... optional clear sky heating rates
1208 
1209  if ( lhsw0 ) then
1210  fnet(1) = flxd0(1) - flxu0(1)
1211 
1212  do k = 2, nlp1
1213  fnet(k) = flxd0(k) - flxu0(k)
1214  hsw0(j1,k-1) = (fnet(k)-fnet(k-1)) * rfdelp(k-1)
1215  enddo
1216  endif
1217 
1218 !! --- ... optional spectral band heating rates
1219 
1220  if ( lhswb ) then
1221  do mb = 1, nbdsw
1222  fnet(1) = fxdnc(1,mb) - fxupc(1,mb)
1223 
1224  do k = 1, nlay
1225  fnet(k+1) = fxdnc(k+1,mb) - fxupc(k+1,mb)
1226  hswb(j1,k,mb) = (fnet(k+1) - fnet(k)) * rfdelp(k)
1227  enddo
1228  enddo
1229  endif
1230 
1231  endif ! if_ivflip
1232 
1233  enddo lab_do_ipt
1234 
1235  return
1236 !...................................

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine module_radsw_main::taumol ( )
private

Definition at line 3593 of file radsw_main.f.

References f_one, module_radsw_sflux::ibx, module_radsw_sflux::ix1, module_radsw_sflux::ix2, module_radsw_sflux::layreffr, module_radsw_parameters::nbhgh, module_radsw_parameters::nblow, module_radsw_parameters::ng, module_radsw_parameters::ngs, nspa, nspb, oneminus, module_radsw_sflux::scalekur, module_radsw_sflux::sfluxref01, module_radsw_sflux::sfluxref02, module_radsw_sflux::sfluxref03, module_radsw_sflux::specwt, module_radsw_sflux::strrat, taumol16(), taumol17(), taumol18(), taumol19(), taumol20(), taumol21(), taumol22(), taumol23(), taumol24(), taumol25(), taumol26(), taumol27(), taumol28(), and taumol29().

Referenced by swrad().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

real (kind=kind_phys), parameter module_radsw_main::amdo3 = con_amd/con_amo3
private

Definition at line 304 of file radsw_main.f.

Referenced by swrad().

304  real (kind=kind_phys), parameter :: amdo3 = con_amd/con_amo3
real (kind=kind_phys), parameter module_radsw_main::amdw = con_amd/con_amw
private

Definition at line 303 of file radsw_main.f.

Referenced by swrad().

303  real (kind=kind_phys), parameter :: amdw = con_amd/con_amw
real (kind=kind_phys), parameter module_radsw_main::bpade = 1.0/0.278
private

Definition at line 294 of file radsw_main.f.

Referenced by rswinit(), spcvrtc(), and spcvrtm().

294  real (kind=kind_phys), parameter :: bpade = 1.0/0.278 ! pade approx constant
real (kind=kind_phys), parameter module_radsw_main::eps = 1.0e-6
private

Definition at line 292 of file radsw_main.f.

Referenced by spcvrtc(), and spcvrtm().

292  real (kind=kind_phys), parameter :: eps = 1.0e-6
real (kind=kind_phys), dimension(0:ntbmx) module_radsw_main::exp_tbl
private

Definition at line 340 of file radsw_main.f.

Referenced by rswinit(), spcvrtc(), and spcvrtm().

340  real (kind=kind_phys) :: exp_tbl(0:ntbmx)
real (kind=kind_phys), parameter module_radsw_main::f_one = 1.0
private

Definition at line 300 of file radsw_main.f.

Referenced by cldprop(), mcica_subcol(), setcoef(), spcvrtc(), spcvrtm(), swflux(), swrad(), taumol(), taumol16(), taumol17(), taumol18(), taumol19(), taumol21(), taumol22(), taumol24(), and taumol28().

300  real (kind=kind_phys), parameter :: f_one = 1.0
real (kind=kind_phys), parameter module_radsw_main::f_zero = 0.0
private

Definition at line 299 of file radsw_main.f.

Referenced by cldprop(), setcoef(), spcvrtc(), spcvrtm(), swflux(), swrad(), taumol23(), and taumol26().

299  real (kind=kind_phys), parameter :: f_zero = 0.0
real (kind=kind_phys), parameter module_radsw_main::ftiny = 1.0e-12
private

Definition at line 296 of file radsw_main.f.

Referenced by cldprop(), spcvrtc(), spcvrtm(), and swrad().

296  real (kind=kind_phys), parameter :: ftiny = 1.0e-12
real (kind=kind_phys) module_radsw_main::heatfac
private

Definition at line 345 of file radsw_main.f.

Referenced by rswinit(), and swrad().

345  real (kind=kind_phys) :: heatfac
integer, dimension(nblow:nbhgh) module_radsw_main::idxebc
private

Definition at line 307 of file radsw_main.f.

Referenced by cldprop().

integer, dimension(nblow:nbhgh) module_radsw_main::idxsfc
private

Definition at line 307 of file radsw_main.f.

Referenced by spcvrtc(), and spcvrtm().

integer, parameter module_radsw_main::ipsdsw0 = 1
private

Definition at line 349 of file radsw_main.f.

Referenced by swrad().

349  integer, parameter :: ipsdsw0 = 1 ! initial permutation seed
logical module_radsw_main::lfdncmp = .false.
private

Definition at line 336 of file radsw_main.f.

Referenced by swrad().

336  logical :: lfdncmp= .false.
logical module_radsw_main::lflxprf = .false.
private

Definition at line 335 of file radsw_main.f.

Referenced by swrad().

335  logical :: lflxprf= .false.
logical module_radsw_main::lhsw0 = .false.
private

Definition at line 334 of file radsw_main.f.

Referenced by swrad().

334  logical :: lhsw0 = .false.
logical module_radsw_main::lhswb = .false.
private

Definition at line 333 of file radsw_main.f.

Referenced by swrad().

333  logical :: lhswb = .false.
integer, dimension(nblow:nbhgh) module_radsw_main::nspa
private

Definition at line 307 of file radsw_main.f.

Referenced by taumol().

307  integer, dimension(nblow:nbhgh) :: nspa, nspb, idxebc, idxsfc
integer, dimension(nblow:nbhgh) module_radsw_main::nspb
private

Definition at line 307 of file radsw_main.f.

Referenced by taumol().

integer, parameter module_radsw_main::nuvb = 27
private

Definition at line 329 of file radsw_main.f.

Referenced by spcvrtc(), and spcvrtm().

329  integer, parameter :: nuvb = 27 !uv-b band index
real (kind=kind_phys), parameter module_radsw_main::oneminus = 1.0 - eps
private

Definition at line 293 of file radsw_main.f.

Referenced by spcvrtc(), spcvrtm(), swrad(), taumol(), taumol16(), taumol17(), taumol18(), taumol19(), taumol21(), taumol22(), taumol24(), and taumol28().

293  real (kind=kind_phys), parameter :: oneminus= 1.0 - eps
real (kind=kind_phys), parameter module_radsw_main::s0 = 1368.22
private

Definition at line 297 of file radsw_main.f.

Referenced by swrad().

297  real (kind=kind_phys), parameter :: s0 = 1368.22 ! internal solar const
real (kind=kind_phys), parameter module_radsw_main::stpfac = 296.0/1013.0
private

Definition at line 295 of file radsw_main.f.

Referenced by setcoef().

295  real (kind=kind_phys), parameter :: stpfac = 296.0/1013.0
character(40), parameter module_radsw_main::vtagsw ='NCEP SW v5.1 Nov 2012 -RRTMG-SW v3.8 '
private

Definition at line 280 of file radsw_main.f.

Referenced by rswinit().

280  character(40), parameter :: &
281  & VTAGSW='NCEP SW v5.1 Nov 2012 -RRTMG-SW v3.8 '