LCOV - code coverage report
Current view: | top level - physics/camrt - radsw.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 606 | 646 | 93.8 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 3 | 3 | 100.0 %

              Line data    Source code
           1             :
           2             : module radsw
           3             : !-----------------------------------------------------------------------
           4             : !
           5             : ! Purpose: Solar radiation calculations.
           6             : !
           7             : !-----------------------------------------------------------------------
           8             : use shr_kind_mod,    only: r8 => shr_kind_r8
           9             : use ppgrid,          only: pcols, pver, pverp
          10             : use cam_abortutils,  only: endrun
          11             : use cam_history,     only: outfld
          12             : use scamMod,         only: single_column,scm_crm_mode,have_asdir, &
          13             :                            asdirobs, have_asdif, asdifobs, have_aldir, &
          14             :                            aldirobs, have_aldif, aldifobs
          15             : use cam_logfile,     only: iulog
          16             : use radconstants,    only: nswbands, get_sw_spectral_boundaries, &
          17             :                            idx_sw_diag, indxsl
          18             :
          19             : implicit none
          20             :
          21             : private
          22             : save
          23             :
          24             : ! Public methods
          25             :
          26             : public ::&
          27             :    radsw_init,   &! initialize constants
          28             :    radcswmx       ! driver for solar radiation code
          29             :
          30             : ! Private module data
          31             :
          32             : real(r8) :: gravit     ! Acceleration of gravity
          33             : real(r8) :: rga        ! 1./gravit
          34             : real(r8) :: sslp       ! Standard sea-level pressure
          35             :
          36             : !===============================================================================
          37             : CONTAINS
          38             : !===============================================================================
          39             :
          40      749232 : subroutine radcswmx(lchnk   ,ncol    ,                         &
          41             :                     E_pint    ,E_pmid    ,E_h2ommr  ,E_o3mmr   , &
          42             :                     E_o2mmr   ,E_cld     ,E_cicewp  ,E_cliqwp  ,E_rel     , &
          43             :                     E_rei     ,eccf      ,E_coszrs  ,solin     , &
          44             :                     E_asdir   ,E_asdif   ,E_aldir   ,E_aldif   ,nmxrgn  , &
          45             :                     pmxrgn  ,qrs,qrsc,fsnt    ,fsntc  ,fsdtoa,  fsntoa,   &
          46             :                     fsutoa ,fsntoac, fsnirtoa,fsnrtoac,fsnrtoaq,fsns    , &
          47             :                     fsnsc   ,fsdsc   ,fsds    ,sols    ,soll    , &
          48             :                     solsd   ,solld   , fns     ,fcns            , &
          49             :                     Nday    ,Nnite   ,IdxDay  ,IdxNite, E_co2mmr, &
          50             :                     E_aer_tau, E_aer_tau_w, E_aer_tau_w_g, E_aer_tau_w_f, tauxcl_out, tauxci_out)
          51             : !-----------------------------------------------------------------------
          52             : !
          53             :
          54             : ! Purpose:
          55             : ! Solar radiation code
          56             : !
          57             : ! Method:
          58             : ! Basic method is Delta-Eddington as described in:
          59             : !
          60             : ! Briegleb, Bruce P., 1992: Delta-Eddington
          61             : ! Approximation for Solar Radiation in the NCAR Community Climate Model,
          62             : ! Journal of Geophysical Research, Vol 97, D7, pp7603-7612).
          63             : !
          64             : ! Five changes to the basic method described above are:
          65             : ! (1) addition of sulfate aerosols (Kiehl and Briegleb, 1993)
          66             : ! (2) the distinction between liquid and ice particle clouds
          67             : ! (Kiehl et al, 1996);
          68             : ! (3) provision for calculating TOA fluxes with spectral response to
          69             : ! match Nimbus-7 visible/near-IR radiometers (Collins, 1998);
          70             : ! (4) max-random overlap (Collins, 2001)
          71             : ! (5) The near-IR absorption by H2O was updated in 2003 by Collins,
          72             : !     Lee-Taylor, and Edwards for consistency with the new line data in
          73             : !     Hitran 2000 and the H2O continuum version CKD 2.4.  Modifications
          74             : !     were optimized by reducing RMS errors in heating rates relative
          75             : !     to a series of benchmark calculations for the 5 standard AFGL
          76             : !     atmospheres.  The benchmarks were performed using DISORT2 combined
          77             : !     with GENLN3.  The near-IR scattering optical depths for Rayleigh
          78             : !     scattering were also adjusted, as well as the correction for
          79             : !     stratospheric heating by H2O.
          80             : !
          81             : ! The treatment of maximum-random overlap is described in the
          82             : ! comment block "INDEX CALCULATIONS FOR MAX OVERLAP".
          83             : !
          84             : ! Divides solar spectrum into 19 intervals from 0.2-5.0 micro-meters.
          85             : ! solar flux fractions specified for each interval. allows for
          86             : ! seasonally and diurnally varying solar input.  Includes molecular,
          87             : ! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud,
          88             : ! and surface absorption. Computes delta-eddington reflections and
          89             : ! transmissions assuming homogeneously mixed layers. Adds the layers
          90             : ! assuming scattering between layers to be isotropic, and distinguishes
          91             : ! direct solar beam from scattered radiation.
          92             : !
          93             : ! Longitude loops are broken into 1 or 2 sections, so that only daylight
          94             : ! (i.e. coszrs > 0) computations are done.
          95             : !
          96             : ! Note that an extra layer above the model top layer is added.
          97             : !
          98             : ! cgs units are used.
          99             : !
         100             : ! Special diagnostic calculation of the clear sky surface and total column
         101             : ! absorbed flux is also done for cloud forcing diagnostics.
         102             : !
         103             : !-----------------------------------------------------------------------
         104             : !
         105             : ! D. Parks (NEC) 09/11/03
         106             : ! Restructuring of routine to support SX vector architecture.
         107             : !
         108             : ! Possible improvements:
         109             : !
         110             : ! 1. Look at vectorizing index calculations for maximum overlap.
         111             : !
         112             : ! 2. Consider making innermost loop in flux computations the number
         113             : !    of spectral intervals.  Given that NS is fixed at 19, the trade-off
         114             : !    will be stride one memory accesses of length 19 versus indirect
         115             : !    addressing (list vector - gather/scatter) with potential vector
         116             : !    lenghts of the number of day light points.  Vectorizing on the number
         117             : !    of spectral intervals seems worthwhile for low resolution models (T42),
         118             : !    but might be inefficient with higher resolutions.
         119             : !
         120             : ! 3. Move the linearization of daylight points (compression/expansion) out
         121             : !    of radcswmx and into d_p_coupling.  This would eliminate the cost of
         122             : !    routines CmpDayNite and ExpDayNite.
         123             : !
         124             : ! 4. Look at expliciting computing all streams in upward propagation of
         125             : !    radiation. There would be additional floating point operations in
         126             : !    exchange for the elimination of indirect addressing.
         127             : !
         128             : !-----------------------------------------------------------------------
         129             :
         130             :    use rad_solar_var,    only: get_variability
         131             :    use cmparray_mod,     only: CmpDayNite, ExpDayNite
         132             :    use quicksort,        only: quick_sort
         133             :    use phys_control,     only: phys_getopts
         134             :    use solar_irrad_data, only: sol_tsi, do_spctrl_scaling, ref_tsi
         135             :    use radconstants,     only: frcsol, ph2o, pco2, po2
         136             :
         137             : !-----------------------Constants for new band (640-700 nm)-------------
         138             :    real(r8) v_raytau_35
         139             :    real(r8) v_raytau_64
         140             :    real(r8) v_abo3_35
         141             :    real(r8) v_abo3_64
         142             :    parameter( &
         143             :         v_raytau_35 = 0.155208_r8, &
         144             :         v_raytau_64 = 0.0392_r8, &
         145             :         v_abo3_35 = 2.4058030e+01_r8, &
         146             :         v_abo3_64 = 2.210e+01_r8 &
         147             :         )
         148             :
         149             :
         150             : !-------------Parameters for accelerating max-random solution-------------
         151             : !
         152             : ! The solution time scales like prod(j:1->N) (1 + n_j) where
         153             : ! N   = number of max-overlap regions (nmxrgn)
         154             : ! n_j = number of unique cloud amounts in region j
         155             : !
         156             : ! Therefore the solution cost can be reduced by decreasing n_j.
         157             : ! cldmin reduces n_j by treating cloud amounts < cldmin as clear sky.
         158             : ! cldeps reduces n_j by treating cloud amounts identical to log(1/cldeps)
         159             : ! decimal places as identical
         160             : !
         161             : ! areamin reduces the cost by dropping configurations that occupy
         162             : ! a surface area < areamin of the model grid box.  The surface area
         163             : ! for a configuration C(j,k_j), where j is the region number and k_j is the
         164             : ! index for a unique cloud amount (in descending order from biggest to
         165             : ! smallest clouds) in region j, is
         166             : !
         167             : ! A = prod(j:1->N) [C(j,k_j) - C(j,k_j+1)]
         168             : !
         169             : ! where C(j,0) = 1.0 and C(j,n_j+1) = 0.0.
         170             : !
         171             : ! nconfgmax reduces the cost and improves load balancing by setting an upper
         172             : ! bound on the number of cloud configurations in the solution.  If the number
         173             : ! of configurations exceeds nconfgmax, the nconfgmax configurations with the
         174             : ! largest area are retained, and the fluxes are normalized by the total area
         175             : ! of these nconfgmax configurations.  For the current max/random overlap
         176             : ! assumption (see subroutine cldovrlap), 30 levels, and cloud-amount
         177             : ! parameterization, the mean and RMS number of configurations are
         178             : ! both roughly 5.  nconfgmax has been set to the mean+2*RMS number, or 15.
         179             : !
         180             : ! Minimum cloud amount (as a fraction of the grid-box area) to
         181             : ! distinguish from clear sky
         182             : !
         183             :    real(r8) cldmin
         184             :    parameter (cldmin = 1.0e-80_r8)
         185             : !
         186             : ! Minimimum horizontal area (as a fraction of the grid-box area) to retain
         187             : ! for a unique cloud configuration in the max-random solution
         188             : !
         189             :    real(r8) areamin
         190             :    parameter (areamin = 0.01_r8)
         191             : !
         192             : ! Decimal precision of cloud amount (0 -> preserve full resolution;
         193             : ! 10^-n -> preserve n digits of cloud amount)
         194             : !
         195             :    real(r8) cldeps
         196             :    parameter (cldeps = 0.0_r8)
         197             : !
         198             : ! Maximum number of configurations to include in solution
         199             : !
         200             :    integer nconfgmax
         201             :    parameter (nconfgmax = 15)
         202             : !------------------------------Commons----------------------------------
         203             : !
         204             : ! Input arguments
         205             : !
         206             :    integer, intent(in) :: lchnk             ! chunk identifier
         207             :    integer, intent(in) :: ncol              ! number of atmospheric columns
         208             :
         209             :    integer,intent(in) :: Nday                      ! Number of daylight columns
         210             :    integer,intent(in) :: Nnite                     ! Number of night columns
         211             :    integer,intent(in), dimension(pcols) :: IdxDay  ! Indicies of daylight coumns
         212             :    integer,intent(in), dimension(pcols) :: IdxNite ! Indicies of night coumns
         213             :
         214             :
         215             :    real(r8), intent(in) :: E_pmid(pcols,pver) ! Level pressure
         216             :    real(r8), intent(in) :: E_pint(pcols,pverp) ! Interface pressure
         217             :    real(r8), intent(in) :: E_h2ommr(pcols,pver) ! Specific humidity (h2o mass mix ratio)
         218             :    real(r8), intent(in) :: E_o3mmr(pcols,pver) ! Ozone mass mixing ratio
         219             :    real(r8), intent(in) :: E_o2mmr(pcols) ! oxygen mass mixing ratio
         220             : !
         221             :    real(r8), intent(in) :: E_cld(pcols,pver)  ! Fractional cloud cover
         222             :    real(r8), intent(in) :: E_cicewp(pcols,pver) ! in-cloud cloud ice water path
         223             :    real(r8), intent(in) :: E_cliqwp(pcols,pver) ! in-cloud cloud liquid water path
         224             :    real(r8), intent(in) :: E_rel(pcols,pver)  ! Liquid effective drop size (microns)
         225             :    real(r8), intent(in) :: E_rei(pcols,pver)  ! Ice effective drop size (microns)
         226             : !
         227             :    real(r8), intent(in) :: eccf             ! Eccentricity factor (1./earth-sun dist^2)
         228             :    real(r8), intent(in) :: E_coszrs(pcols)    ! Cosine solar zenith angle
         229             :    real(r8), intent(in) :: E_asdir(pcols)     ! 0.2-0.7 micro-meter srfc alb: direct rad
         230             :    real(r8), intent(in) :: E_aldir(pcols)     ! 0.7-5.0 micro-meter srfc alb: direct rad
         231             :    real(r8), intent(in) :: E_asdif(pcols)     ! 0.2-0.7 micro-meter srfc alb: diffuse rad
         232             :    real(r8), intent(in) :: E_aldif(pcols)     ! 0.7-5.0 micro-meter srfc alb: diffuse rad
         233             :    real(r8), intent(in) :: E_co2mmr(pcols)    ! co2 column mean mmr
         234             :
         235             : !
         236             : ! Aerosol radiative property arrays
         237             : !
         238             :    real(r8),intent(in) :: E_aer_tau    (pcols,0:pver,nswbands) ! aerosol extinction optical depth
         239             :    real(r8),intent(in) :: E_aer_tau_w  (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
         240             :    real(r8),intent(in) :: E_aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
         241             :    real(r8),intent(in) :: E_aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
         242             :
         243             : !
         244             : ! IN/OUT arguments
         245             : !
         246             :    real(r8), intent(inout) :: pmxrgn(pcols,pverp) ! Maximum values of pressure for each
         247             : !                                                 !    maximally overlapped region.
         248             : !                                                 !    0->pmxrgn(i,1) is range of pressure for
         249             : !                                                 !    1st region,pmxrgn(i,1)->pmxrgn(i,2) for
         250             : !                                                 !    2nd region, etc
         251             :    integer, intent(inout) ::  nmxrgn(pcols)    ! Number of maximally overlapped regions
         252             : !
         253             : ! Output arguments
         254             : !
         255             :
         256             :    real(r8), intent(out) :: solin(pcols)     ! Incident solar flux
         257             :    real(r8), intent(out) :: qrs (pcols,pver) ! Solar heating rate
         258             :    real(r8), intent(out) :: qrsc(pcols,pver) ! Clearsky solar heating rate
         259             :    real(r8), intent(out) :: fsns(pcols)      ! Surface absorbed solar flux
         260             :    real(r8), intent(out) :: fsnt(pcols)      ! Total column absorbed solar flux
         261             :    real(r8), intent(out) :: fsntoa(pcols)    ! Net solar flux at TOA
         262             :    real(r8), intent(out) :: fsutoa(pcols)    ! Upward solar flux at TOA
         263             :    real(r8), intent(out) :: fsds(pcols)      ! Flux shortwave downwelling surface
         264             : !
         265             :    real(r8), intent(out) :: fsnsc(pcols)     ! Clear sky surface absorbed solar flux
         266             :    real(r8), intent(out) :: fsdsc(pcols)     ! Clear sky surface downwelling solar flux
         267             :    real(r8), intent(out) :: fsntc(pcols)     ! Clear sky total column absorbed solar flx
         268             :    real(r8), intent(out) :: fsdtoa(pcols)    ! Downwelling solar flux at TOA
         269             :    real(r8), intent(out) :: fsntoac(pcols)   ! Clear sky net solar flx at TOA
         270             :    real(r8), intent(out) :: sols(pcols)      ! Direct solar rad on surface (< 0.7)
         271             :    real(r8), intent(out) :: soll(pcols)      ! Direct solar rad on surface (>= 0.7)
         272             :    real(r8), intent(out) :: solsd(pcols)     ! Diffuse solar rad on surface (< 0.7)
         273             :    real(r8), intent(out) :: solld(pcols)     ! Diffuse solar rad on surface (>= 0.7)
         274             :    real(r8), intent(out) :: fsnirtoa(pcols)  ! Near-IR flux absorbed at toa
         275             :    real(r8), intent(out) :: fsnrtoac(pcols)  ! Clear sky near-IR flux absorbed at toa
         276             :    real(r8), intent(out) :: fsnrtoaq(pcols)  ! Net near-IR flux at toa >= 0.7 microns
         277             :
         278             :    real(r8), intent(out) :: fns(pcols,pverp)   ! net flux at interfaces
         279             :    real(r8), intent(out) :: fcns(pcols,pverp)  ! net clear-sky flux at interfaces
         280             :
         281             :    real(r8), intent(out) :: tauxcl_out(pcols,pver) ! liquid cloud visible sw optical depth
         282             :    real(r8), intent(out) :: tauxci_out(pcols,pver) ! ice cloud visible sw optical depth
         283             :
         284             : !
         285             : !---------------------------Local variables-----------------------------
         286             : !
         287             : ! Local and reordered copies of the intent(in) variables
         288             : !
         289             :    real(r8):: aer_tau    (pcols,0:pver,nswbands) ! aerosol extinction optical depth
         290             :    real(r8):: aer_tau_w  (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
         291             :    real(r8):: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
         292             :    real(r8):: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
         293             :    real(r8) :: pmid(pcols,pver) ! Level pressure
         294             :    real(r8) :: pint(pcols,pverp) ! Interface pressure
         295             :    real(r8) :: h2ommr(pcols,pver) ! Specific humidity (h2o mass mix ratio)
         296             :    real(r8) :: o3mmr(pcols,pver) ! Ozone mass mixing ratio
         297             : !
         298             :    real(r8) :: cld(pcols,pver)  ! Fractional cloud cover
         299             :    real(r8) :: cicewp(pcols,pver) ! in-cloud cloud ice water path
         300             :    real(r8) :: cliqwp(pcols,pver) ! in-cloud cloud liquid water path
         301             :    real(r8) :: rel(pcols,pver)  ! Liquid effective drop size (microns)
         302             :    real(r8) :: rei(pcols,pver)  ! Ice effective drop size (microns)
         303             : !
         304             :    real(r8) :: coszrs(pcols)    ! Cosine solar zenith angle
         305             :    real(r8) :: asdir(pcols)     ! 0.2-0.7 micro-meter srfc alb: direct rad
         306             :    real(r8) :: aldir(pcols)     ! 0.7-5.0 micro-meter srfc alb: direct rad
         307             :    real(r8) :: asdif(pcols)     ! 0.2-0.7 micro-meter srfc alb: diffuse rad
         308             :    real(r8) :: aldif(pcols)     ! 0.7-5.0 micro-meter srfc alb: diffuse rad
         309             :    real(r8) :: co2mmr(pcols)    ! co2 column mean mmr
         310             :    real(r8) :: o2mmr(pcols)     ! o2 column mean mmr
         311             :
         312             :    real(r8) :: tot_irrad
         313             : !
         314             : ! Max/random overlap variables
         315             : !
         316             :    real(r8) asort(pverp)     ! 1 - cloud amounts to be sorted for max ovrlp.
         317             :    real(r8) atmp             ! Temporary storage for sort when nxs = 2
         318             :    real(r8) cld0             ! 1 - (cld amt) used to make wstr, cstr, nstr
         319             :    real(r8) totwgt(pcols)    ! Total of xwgts = total fractional area of
         320             : !   grid-box covered by cloud configurations
         321             : !   included in solution to fluxes
         322             :
         323             :    real(r8) wgtv(nconfgmax)  ! Weights for fluxes
         324             : !   1st index is configuration number
         325             :    real(r8) wstr(pverp,pverp) ! area weighting factors for streams
         326             : !   1st index is for stream #,
         327             : !   2nd index is for region #
         328             :
         329             :    real(r8) xexpt            ! solar direct beam trans. for layer above
         330             :    real(r8) xrdnd            ! diffuse reflectivity for layer above
         331             :    real(r8) xrupd            ! diffuse reflectivity for layer below
         332             :    real(r8) xrups            ! direct-beam reflectivity for layer below
         333             :    real(r8) xtdnt            ! total trans for layers above
         334             :
         335             :    real(r8) xwgt             ! product of cloud amounts
         336             :
         337             :    real(r8) yexpt            ! solar direct beam trans. for layer above
         338             :    real(r8) yrdnd            ! diffuse reflectivity for layer above
         339             :    real(r8) yrupd            ! diffuse reflectivity for layer below
         340             :    real(r8) ytdnd            ! dif-beam transmission for layers above
         341             :    real(r8) ytupd            ! dif-beam transmission for layers below
         342             :
         343             :    real(r8) zexpt            ! solar direct beam trans. for layer above
         344             :    real(r8) zrdnd            ! diffuse reflectivity for layer above
         345             :    real(r8) zrupd            ! diffuse reflectivity for layer below
         346             :    real(r8) zrups            ! direct-beam reflectivity for layer below
         347             :    real(r8) ztdnt            ! total trans for layers above
         348             :
         349             :    logical new_term          ! Flag for configurations to include in fluxes
         350             :    logical region_found      ! flag for identifying regions
         351             :
         352             :    integer ccon(nconfgmax,0:pverp,pcols)
         353             : ! flags for presence of clouds
         354             : !   1st index is for level # (including
         355             : !    layer above top of model and at surface)
         356             : !   2nd index is for configuration #
         357             :    integer cstr(0:pverp,pverp)
         358             : ! flags for presence of clouds
         359             : !   1st index is for level # (including
         360             : !    layer above top of model and at surface)
         361             : !   2nd index is for stream #
         362             :    integer icond(nconfgmax,0:pverp,pcols)
         363             : ! Indices for copying rad. properties from
         364             : !     one identical downward cld config.
         365             : !     to another in adding method (step 2)
         366             : !   1st index is for interface # (including
         367             : !     layer above top of model and at surface)
         368             : !   2nd index is for configuration # range
         369             :    integer iconu(nconfgmax,0:pverp,pcols)
         370             : ! Indices for copying rad. properties from
         371             : !     one identical upward configuration
         372             : !     to another in adding method (step 2)
         373             : !   1st index is for interface # (including
         374             : !     layer above top of model and at surface)
         375             : !   2nd index is for configuration # range
         376             :    integer iconfig           ! Counter for random-ovrlap configurations
         377             :    integer irgn              ! Index for max-overlap regions
         378             :    integer is0               ! Lower end of stream index range
         379             :    integer is1               ! Upper end of stream index range
         380             :    integer isn               ! Stream index
         381             :    integer istr(pverp+1)     ! index for stream #s during flux calculation
         382             :    integer istrtd(0:nconfgmax+1,0:pverp,pcols)
         383             : ! indices into icond
         384             : !   1st index is for interface # (including
         385             : !     layer above top of model and at surface)
         386             : !   2nd index is for configuration # range
         387             :    integer istrtu(0:nconfgmax+1,0:pverp,pcols)
         388             : ! indices into iconu
         389             : !   1st index is for interface # (including
         390             : !     layer above top of model and at surface)
         391             : !   2nd index is for configuration # range
         392             :    integer j                 ! Configuration index
         393             :    integer jj                ! Configuration index
         394             :    integer k1                ! Level index
         395             :    integer k2                ! Level index
         396             :    integer ksort(pverp)      ! Level indices of cloud amounts to be sorted
         397             :    integer ktmp              ! Temporary storage for sort when nxs = 2
         398             :    integer kx1(0:pverp)      ! Level index for top of max-overlap region
         399             :    integer kx2(0:pverp)      ! Level index for bottom of max-overlap region
         400             :    integer l                 ! Index
         401             :    integer l0                ! Index
         402             :    integer mrgn              ! Counter for nrgn
         403             :    integer mstr              ! Counter for nstr
         404             :    integer n0                ! Number of configurations with ccon(:,k,:)==0
         405             :    integer n1                ! Number of configurations with ccon(:,k,:)==1
         406             :    integer nconfig(pcols)    ! Number of random-ovrlap configurations
         407             :    integer nconfigm          ! Value of config before testing for areamin,
         408             : !    nconfgmax
         409             :    integer npasses           ! number of passes over the indexing loop
         410             :    integer nrgn              ! Number of max overlap regions at current
         411             : !    longitude
         412             :    integer nstr(pverp)       ! Number of unique cloud configurations
         413             : !   ("streams") in a max-overlapped region
         414             : !   1st index is for region #
         415             :    integer nuniq             ! # of unique cloud configurations
         416             :    integer nuniqd(0:pverp,pcols)   ! # of unique cloud configurations: TOA
         417             : !   to level k
         418             :    integer nuniqu(0:pverp,pcols)   ! # of unique cloud configurations: surface
         419             : !   to level k
         420             :    integer nxs               ! Number of cloudy layers between k1 and k2
         421             :    integer ptr0(nconfgmax)   ! Indices of configurations with ccon(:,k,:)==0
         422             :    integer ptr1(nconfgmax)   ! Indices of configurations with ccon(:,k,:)==1
         423             :    integer ptrc(nconfgmax)   ! Pointer for configurations sorted by wgtv
         424             :    integer, dimension(1) :: min_idx  ! required for return val of func minloc
         425             :
         426             : !
         427             : ! Other
         428             : !
         429             :    integer ns                ! Spectral loop index
         430             :    integer i                 ! Longitude loop index
         431             :    integer k                 ! Level loop index
         432             :    integer km1               ! k - 1
         433             :    integer kp1               ! k + 1
         434             :    integer n                 ! Loop index for daylight
         435             :    integer ksz               ! dust size bin index
         436             :    integer kaer              ! aerosol group index
         437             : !
         438             : ! A. Slingo's data for cloud particle radiative properties (from 'A GCM
         439             : ! Parameterization for the Shortwave Properties of Water Clouds' JAS
         440             : ! vol. 46 may 1989 pp 1419-1427)
         441             : !
         442             :
         443             :    real(r8), parameter :: abarl(4) = (/ 2.817e-02_r8, 2.682e-02_r8,2.264e-02_r8,1.281e-02_r8/)
         444             :    real(r8), parameter :: bbarl(4) = (/ 1.305_r8    , 1.346_r8    ,1.454_r8    ,1.641_r8    /)
         445             :    real(r8), parameter :: cbarl(4) = (/-5.62e-08_r8 ,-6.94e-06_r8 ,4.64e-04_r8 ,0.201_r8    /)
         446             :    real(r8), parameter :: dbarl(4) = (/ 1.63e-07_r8 , 2.35e-05_r8 ,1.24e-03_r8 ,7.56e-03_r8 /)
         447             :    real(r8), parameter :: ebarl(4) = (/ 0.829_r8    , 0.794_r8    ,0.754_r8    ,0.826_r8    /)
         448             :    real(r8), parameter :: fbarl(4) = (/ 2.482e-03_r8, 4.226e-03_r8,6.560e-03_r8,4.353e-03_r8/)
         449             :
         450             :    real(r8) :: abarli           ! A coefficient for current spectral band
         451             :    real(r8) :: bbarli           ! B coefficient for current spectral band
         452             :    real(r8) :: cbarli           ! C coefficient for current spectral band
         453             :    real(r8) :: dbarli           ! D coefficient for current spectral band
         454             :    real(r8) :: ebarli           ! E coefficient for current spectral band
         455             :    real(r8) :: fbarli           ! F coefficient for current spectral band
         456             : !
         457             : ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor
         458             : ! greater than 20 micro-meters
         459             : !
         460             : ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
         461             : !
         462             :    real(r8), parameter :: abari(4) = (/ 3.448e-03_r8, 3.448e-03_r8,3.448e-03_r8,3.448e-03_r8/)
         463             :    real(r8), parameter :: bbari(4) = (/ 2.431_r8    , 2.431_r8    ,2.431_r8    ,2.431_r8    /)
         464             :    real(r8), parameter :: cbari(4) = (/ 1.00e-05_r8 , 1.10e-04_r8 ,1.861e-02_r8,.46658_r8   /)
         465             :    real(r8), parameter :: dbari(4) = (/ 0.0_r8      , 1.405e-05_r8,8.328e-04_r8,2.05e-05_r8 /)
         466             :    real(r8), parameter :: ebari(4) = (/ 0.7661_r8   , 0.7730_r8   ,0.794_r8    ,0.9595_r8   /)
         467             :    real(r8), parameter :: fbari(4) = (/ 5.851e-04_r8, 5.665e-04_r8,7.267e-04_r8,1.076e-04_r8/)
         468             :
         469             :    real(r8) :: abarii           ! A coefficient for current spectral band
         470             :    real(r8) :: bbarii           ! B coefficient for current spectral band
         471             :    real(r8) :: cbarii           ! C coefficient for current spectral band
         472             :    real(r8) :: dbarii           ! D coefficient for current spectral band
         473             :    real(r8) :: ebarii           ! E coefficient for current spectral band
         474             :    real(r8) :: fbarii           ! F coefficient for current spectral band
         475             :
         476             : !
         477             : ! UPDATE TO H2O NEAR-IR: Delta optimized for Hitran 2K and CKD 2.4
         478             : !
         479             :    real(r8), parameter :: delta = 0.0014257179260883_r8
         480             : !
         481             : ! END UPDATE
         482             : !
         483             :    real(r8) :: albdir(pcols,nswbands) ! Current spc intrvl srf alb to direct rad
         484             :    real(r8) :: albdif(pcols,nswbands) ! Current spc intrvl srf alb to diffuse rad
         485             : !
         486             : ! Next series depends on spectral interval
         487             : !
         488             :    real(r8) :: wgtint           ! Weight for specific spectral interval
         489             :
         490             : !
         491             : ! weight for 0.64 - 0.7 microns  appropriate to clear skies over oceans
         492             : !
         493             :    real(r8), parameter :: nirwgt(nswbands) = &
         494             :               (/  0.0_r8,   0.0_r8,   0.0_r8,      0.0_r8,   0.0_r8, &
         495             :                   0.0_r8,   0.0_r8,   0.0_r8, 0.320518_r8,   1.0_r8,  1.0_r8, &
         496             :                   1.0_r8,   1.0_r8,   1.0_r8,      1.0_r8,   1.0_r8, &
         497             :                   1.0_r8,   1.0_r8,   1.0_r8 /)
         498             :
         499             : !
         500             : ! UPDATE TO H2O NEAR-IR: Rayleigh scattering optimized for Hitran 2K & CKD 2.4
         501             : !
         502             :    real(r8), parameter :: raytau(nswbands) = &
         503             :                (/ 4.020_r8, 2.180_r8, 1.700_r8, 1.450_r8, 1.250_r8, &
         504             :                   1.085_r8, 0.730_r8, v_raytau_35, v_raytau_64, &
         505             :                   0.02899756_r8, 0.01356763_r8, 0.00537341_r8, &
         506             :                   0.00228515_r8, 0.00105028_r8, 0.00046631_r8, &
         507             :                   0.00025734_r8, &
         508             :                  .0001_r8, .0001_r8, .0001_r8/)
         509             : !
         510             : ! END UPDATE
         511             : !
         512             :
         513             : !
         514             : ! Absorption coefficients
         515             : !
         516             : !
         517             : ! UPDATE TO H2O NEAR-IR: abh2o optimized for Hitran 2K and CKD 2.4
         518             : !
         519             :    real(r8), parameter :: abh2o(nswbands) = &
         520             :              (/    .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         521             :                    .000_r8,     .000_r8,    .000_r8,    .000_r8,    &
         522             :                    0.00256608_r8,  0.06310504_r8,   0.42287445_r8, 2.45397941_r8, &
         523             :                   11.20070807_r8, 47.66091389_r8, 240.19010243_r8, &
         524             :                    .000_r8,    .000_r8,    .000_r8/)
         525             : !
         526             : ! END UPDATE
         527             : !
         528             :
         529             :    real(r8), parameter :: abo3(nswbands) = &
         530             :              (/5.370e+04_r8, 13.080e+04_r8,  9.292e+04_r8, 4.530e+04_r8, 1.616e+04_r8, &
         531             :                4.441e+03_r8,  1.775e+02_r8, v_abo3_35, v_abo3_64,      .000_r8, &
         532             :                .000_r8,   .000_r8    ,   .000_r8   ,   .000_r8   ,      .000_r8, &
         533             :                .000_r8,   .000_r8    ,   .000_r8   ,   .000_r8    /)
         534             :
         535             :    real(r8), parameter :: abco2(nswbands) = &
         536             :               (/   .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         537             :                    .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         538             :                    .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         539             :                    .000_r8,     .094_r8,    .196_r8,   1.963_r8/)
         540             :
         541             :    real(r8), parameter :: abo2(nswbands) = &
         542             :              (/    .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         543             :                    .000_r8,     .000_r8,    .000_r8,1.11e-05_r8,6.69e-05_r8, &
         544             :                    .000_r8,     .000_r8,    .000_r8,    .000_r8,    .000_r8, &
         545             :                    .000_r8,     .000_r8,    .000_r8,    .000_r8/)
         546             : !
         547             : ! Diagnostic and accumulation arrays; note that fswup, and
         548             : ! fswdn are not used in the computation,but are retained for future use.
         549             : !
         550             :    real(r8) solflx(pcols)    ! Solar flux in current interval
         551             :    real(r8) totfld (pcols,0:pver)  ! Spectrally summed flux divergence
         552             :    real(r8) totfldc(pcols,0:pver)  ! Spectrally summed flux divergence (clearsky)
         553             :    real(r8) fswup(pcols,0:pverp)   ! Spectrally summed up flux
         554             :    real(r8) fswdn(pcols,0:pverp)   ! Spectrally summed down flux
         555             : !
         556             : ! Cloud radiative property arrays
         557             : !
         558             :    real(r8) tauxcl(pcols,0:pver) ! water cloud extinction optical depth
         559             :    real(r8) tauxci(pcols,0:pver) ! ice cloud extinction optical depth
         560             :    real(r8) wcl(pcols,0:pver) ! liquid cloud single scattering albedo
         561             :    real(r8) gcl(pcols,0:pver) ! liquid cloud asymmetry parameter
         562             :    real(r8) fcl(pcols,0:pver) ! liquid cloud forward scattered fraction
         563             :    real(r8) wci(pcols,0:pver) ! ice cloud single scattering albedo
         564             :    real(r8) gci(pcols,0:pver) ! ice cloud asymmetry parameter
         565             :    real(r8) fci(pcols,0:pver) ! ice cloud forward scattered fraction
         566             :
         567             : !
         568             : ! Various arrays and other constants:
         569             : !
         570             :    real(r8) pflx(pcols,0:pverp) ! Interface press, including extra layer
         571             :    real(r8) zenfac(pcols)    ! Square root of cos solar zenith angle
         572             :    real(r8) sqrco2(pcols)    ! Square root of the co2 mass mixg ratio
         573             :    real(r8) tmp1             ! Temporary constant array
         574             :    real(r8) tmp2             ! Temporary constant array
         575             :    real(r8) pdel             ! Pressure difference across layer
         576             :    real(r8) path             ! Mass path of layer
         577             :    real(r8) ptop             ! Lower interface pressure of extra layer
         578             :    real(r8) ptho2            ! Used to compute mass path of o2
         579             :    real(r8) ptho3            ! Used to compute mass path of o3
         580             :    real(r8) pthco2           ! Used to compute mass path of co2
         581             :    real(r8) pthh2o           ! Used to compute mass path of h2o
         582             :    real(r8) h2ostr           ! Inverse sq. root h2o mass mixing ratio
         583             :
         584             :    real(r8) wavmin(nswbands) ! Spectral interval minimum wavelength
         585             :    real(r8) wavmax(nswbands) ! Spectral interval maximum wavelength
         586             :    real(r8) wavmid(nswbands) ! Spectral interval middle wavelength
         587             :    real(r8) trayoslp         ! Rayleigh optical depth/standard pressure
         588             :    real(r8) tmp1l            ! Temporary constant array
         589             :    real(r8) tmp2l            ! Temporary constant array
         590             :    real(r8) tmp3l            ! Temporary constant array
         591             :    real(r8) tmp1i            ! Temporary constant array
         592             :    real(r8) tmp2i            ! Temporary constant array
         593             :    real(r8) tmp3i            ! Temporary constant array
         594             :    real(r8) rdenom           ! Multiple scattering term
         595             :    real(r8) rdirexp          ! layer direct ref times exp transmission
         596             :    real(r8) tdnmexp          ! total transmission - exp transmission
         597             :    real(r8) psf(nswbands)      ! Frac of solar flux in spect interval
         598             : !
         599             : ! Layer absorber amounts; note that 0 refers to the extra layer added
         600             : ! above the top model layer
         601             : !
         602             :    real(r8) uh2o(pcols,0:pver) ! Layer absorber amount of h2o
         603             :    real(r8) uo3(pcols,0:pver) ! Layer absorber amount of  o3
         604             :    real(r8) uco2(pcols,0:pver) ! Layer absorber amount of co2
         605             :    real(r8) uo2(pcols,0:pver) ! Layer absorber amount of  o2
         606             : !
         607             : ! Total column absorber amounts:
         608             : !
         609             :    real(r8) uth2o(pcols)     ! Total column  absorber amount of  h2o
         610             :    real(r8) uto3(pcols)      ! Total column  absorber amount of  o3
         611             :    real(r8) utco2(pcols)     ! Total column  absorber amount of  co2
         612             :    real(r8) uto2(pcols)      ! Total column  absorber amount of  o2
         613             : !
         614             : ! These arrays are defined for pver model layers; 0 refers to the extra
         615             : ! layer on top:
         616             : !
         617             :    real(r8) rdir(nswbands,pcols,0:pver) ! Layer reflectivity to direct rad
         618             :    real(r8) rdif(nswbands,pcols,0:pver) ! Layer reflectivity to diffuse rad
         619             :    real(r8) tdir(nswbands,pcols,0:pver) ! Layer transmission to direct rad
         620             :    real(r8) tdif(nswbands,pcols,0:pver) ! Layer transmission to diffuse rad
         621             :    real(r8) explay(nswbands,pcols,0:pver) ! Solar beam exp trans. for layer
         622             :
         623             :    real(r8) rdirc(nswbands,pcols,0:pver) ! Clear Layer reflec. to direct rad
         624             :    real(r8) rdifc(nswbands,pcols,0:pver) ! Clear Layer reflec. to diffuse rad
         625             :    real(r8) tdirc(nswbands,pcols,0:pver) ! Clear Layer trans. to direct rad
         626             :    real(r8) tdifc(nswbands,pcols,0:pver) ! Clear Layer trans. to diffuse rad
         627             :    real(r8) explayc(nswbands,pcols,0:pver) ! Solar beam exp trans. clear layer
         628             :
         629             :    real(r8) fus(pcols,pverp)   ! Upward flux (added for CRM)
         630             :    real(r8) fds(pcols,pverp)   ! Downward flux (added for CRM)
         631             :    real(r8) fusc(pcols,pverp)  ! Upward clear-sky flux (added for CRM)
         632             :    real(r8) fdsc(pcols,pverp) ! Downward clear-sky flux (added for CRM)
         633             :
         634             :    real(r8) flxdiv           ! Flux divergence for layer
         635             :
         636             : !
         637             : ! Temporary arrays for either clear or cloudy values.
         638             : !
         639             :    real(r8), dimension(nswbands) :: Trdir
         640             :    real(r8), dimension(nswbands) :: Trdif
         641             :    real(r8), dimension(nswbands) :: Ttdir
         642             :    real(r8), dimension(nswbands) :: Ttdif
         643             :    real(r8), dimension(nswbands) :: Texplay
         644             : !
         645             : !
         646             : ! Radiative Properties:
         647             : !
         648             : ! There are 1 classes of properties:
         649             : ! (1. All-sky bulk properties
         650             : ! (2. Clear-sky properties
         651             : !
         652             : ! The first set of properties are generated during step 2 of the solution.
         653             : !
         654             : ! These arrays are defined at model interfaces; in 1st index (for level #),
         655             : ! 0 is the top of the extra layer above the model top, and
         656             : ! pverp is the earth surface.  2nd index is for cloud configuration
         657             : ! defined over a whole column.
         658             : !
         659             :    real(r8) exptdn(nswbands,0:pverp,nconfgmax,pcols) ! Sol. beam trans from layers above
         660             :    real(r8) rdndif(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dif rad for layers above
         661             :    real(r8) rupdif(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dif rad for layers below
         662             :    real(r8) rupdir(nswbands,0:pverp,nconfgmax,pcols) ! Ref to dir rad for layers below
         663             :    real(r8) tdntot(nswbands,0:pverp,nconfgmax,pcols) ! Total trans for layers above
         664             : !
         665             : ! Bulk properties used during the clear-sky calculation.
         666             : !
         667             :    real(r8) exptdnc(pcols,0:pverp) ! clr: Sol. beam trans from layers above
         668             :    real(r8) rdndifc(pcols,0:pverp) ! clr: Ref to dif rad for layers above
         669             :    real(r8) rupdifc(pcols,0:pverp) ! clr: Ref to dif rad for layers below
         670             :    real(r8) rupdirc(pcols,0:pverp) ! clr: Ref to dir rad for layers below
         671             :    real(r8) tdntotc(pcols,0:pverp) ! clr: Total trans for layers above
         672             :
         673             :    real(r8) fluxup(nswbands,0:pverp,pcols)  ! Up   flux at model interface
         674             :    real(r8) fluxdn(nswbands,0:pverp,pcols)  ! Down flux at model interface
         675             :    real(r8) wexptdn(nswbands,pcols)   ! Direct solar beam trans. to surface
         676             : !
         677             : ! Scalars used in vectorization
         678             : !
         679             :   integer :: kk
         680             : !
         681             : ! Arrays used in vectorization
         682             : !
         683             :    real(r8) v_wgtv(nconfgmax,pcols)  ! Weights for fluxes
         684             :
         685             :    real(r8) :: rdiff, ro, rn
         686             :    rdiff(ro,rn) = abs((ro-rn)/merge(ro,1.0_r8,ro /= 0.0_r8))
         687             :
         688             : !  solar variability factor
         689             :    real(r8) :: sfac(nswbands)
         690             :
         691             :    character(len=16) :: microp_scheme  ! microphysics scheme
         692             :
         693             :    !-----------------------------------------------------------------------
         694             :    ! START OF CALCULATION
         695             :    !-----------------------------------------------------------------------
         696             :
         697      749232 :    call phys_getopts(microp_scheme_out=microp_scheme)
         698             :
         699             : !
         700             : ! Initialize output fields:
         701             : !
         702    12510432 :    fsds(1:ncol)     = 0.0_r8
         703             :
         704    12510432 :    fsnirtoa(1:ncol) = 0.0_r8
         705    12510432 :    fsnrtoac(1:ncol) = 0.0_r8
         706    12510432 :    fsnrtoaq(1:ncol) = 0.0_r8
         707             :
         708    12510432 :    fsns(1:ncol)     = 0.0_r8
         709    12510432 :    fsnsc(1:ncol)    = 0.0_r8
         710    12510432 :    fsdsc(1:ncol)    = 0.0_r8
         711             :
         712    12510432 :    fsnt(1:ncol)     = 0.0_r8
         713    12510432 :    fsntc(1:ncol)    = 0.0_r8
         714    12510432 :    fsntoa(1:ncol)   = 0.0_r8
         715    12510432 :    fsdtoa(1:ncol)   = 0.0_r8
         716    12510432 :    fsutoa(1:ncol)   = 0.0_r8
         717    12510432 :    fsntoac(1:ncol)  = 0.0_r8
         718             :
         719    12510432 :    solin(1:ncol)    = 0.0_r8
         720             :
         721    12510432 :    sols(1:ncol)     = 0.0_r8
         722    12510432 :    soll(1:ncol)     = 0.0_r8
         723    12510432 :    solsd(1:ncol)    = 0.0_r8
         724    12510432 :    solld(1:ncol)    = 0.0_r8
         725             :
         726   326020464 :    qrs (1:ncol,1:pver) = 0.0_r8
         727   326020464 :    qrsc(1:ncol,1:pver) = 0.0_r8
         728   338530896 :    fns(1:ncol,1:pverp) = 0.0_r8
         729   338530896 :    fcns(1:ncol,1:pverp) = 0.0_r8
         730      749232 :    if (single_column.and.scm_crm_mode) then
         731           0 :       fus(1:ncol,1:pverp) = 0.0_r8
         732           0 :       fds(1:ncol,1:pverp) = 0.0_r8
         733           0 :       fusc(:ncol,:pverp) = 0.0_r8
         734           0 :       fdsc(:ncol,:pverp) = 0.0_r8
         735             :    endif
         736             :
         737      749232 :    tauxcl_out(1:pcols,1:pver) = 0.0_r8
         738      749232 :    tauxci_out(1:pcols,1:pver) = 0.0_r8
         739             :
         740             : !
         741             : ! If night everywhere, return:
         742             : !
         743      749232 :    if ( Nday == 0 ) then
         744             :      return
         745             :    endif
         746             :
         747             : !
         748             : ! Rearrange input arrays
         749             : !
         750             :
         751      391783 :    call CmpDayNite(E_pmid, pmid,        Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         752      391783 :    call CmpDayNite(E_pint, pint,        Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
         753      391783 :    call CmpDayNite(E_h2ommr, h2ommr,    Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         754      391783 :    call CmpDayNite(E_o3mmr, o3mmr,      Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         755      391783 :    call CmpDayNite(E_cld, cld,          Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         756      391783 :    call CmpDayNite(E_cicewp, cicewp,    Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         757      391783 :    call CmpDayNite(E_cliqwp, cliqwp,    Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         758      391783 :    call CmpDayNite(E_rel, rel,          Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         759      391783 :    call CmpDayNite(E_rei, rei,          Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
         760      391783 :    call CmpDayNite(E_coszrs, coszrs,    Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         761      391783 :    call CmpDayNite(E_asdir, asdir,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         762      391783 :    call CmpDayNite(E_aldir, aldir,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         763      391783 :    call CmpDayNite(E_asdif, asdif,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         764      391783 :    call CmpDayNite(E_aldif, aldif,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         765      391783 :    call CmpDayNite(E_co2mmr, co2mmr,    Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         766      391783 :    call CmpDayNite(E_o2mmr, o2mmr,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         767             :
         768      391783 :    call CmpDayNite(pmxrgn,                     Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
         769      391783 :    call CmpDayNite(nmxrgn,                     Nday, IdxDay, Nnite, IdxNite, 1, pcols)
         770      391783 :    call CmpDayNite(E_aer_tau,     aer_tau,     Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
         771      391783 :    call CmpDayNite(E_aer_tau_w,   aer_tau_w,   Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
         772      391783 :    call CmpDayNite(E_aer_tau_w_g, aer_tau_w_g, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
         773      391783 :    call CmpDayNite(E_aer_tau_w_f, aer_tau_w_f, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 0, pver, 1,nswbands)
         774             :
         775      391783 :    if (scm_crm_mode) then
         776             :    ! overwrite albedos for CRM
         777           0 :       if(have_asdir) asdir = asdirobs(1)
         778           0 :       if(have_asdif) asdif = asdifobs(1)
         779           0 :       if(have_aldir) aldir = aldirobs(1)
         780           0 :       if(have_aldif) aldif = aldifobs(1)
         781             :    endif
         782             : !
         783             : ! Perform other initializations
         784             : !
         785      391783 :    tmp1   = 0.5_r8/(gravit*sslp)
         786      391783 :    tmp2   = delta/gravit
         787             :
         788     6272383 :    do i = 1, Nday
         789     6272383 :       sqrco2(i) = sqrt(co2mmr(i))
         790             :    end do
         791             :
         792      391783 :    if ( do_spctrl_scaling ) then
         793           0 :       call get_variability(sfac)
         794           0 :       tot_irrad = ref_tsi
         795             :    else
         796      391783 :       tot_irrad = sol_tsi
         797             :    endif
         798             :
         799    10969924 :    do k=1,pverp
         800   169746124 :       do i=1,Nday
         801   169354341 :          pflx(i,k) = pint(i,k)
         802             :       end do
         803             :    end do
         804             :
         805     6272383 :    do i=1,Nday
         806             : !
         807             : ! Define solar incident radiation and interface pressures:
         808             : !
         809             :
         810     5880600 :          solin(i)  = tot_irrad*1.e3_r8*eccf*coszrs(i)
         811             :
         812     5880600 :          pflx(i,0) = 0._r8
         813             : !
         814             : ! Compute optical paths:
         815             : !
         816     5880600 :          ptop      = pflx(i,1)
         817     5880600 :          ptho2     = o2mmr(i) * ptop / gravit
         818     5880600 :          ptho3     = o3mmr(i,1) * ptop / gravit
         819     5880600 :          pthco2    = sqrco2(i) * (ptop / gravit)
         820     5880600 :          h2ostr    = sqrt( 1._r8 / h2ommr(i,1) )
         821     5880600 :          zenfac(i) = sqrt(coszrs(i))
         822             :          pthh2o    = ptop**2*tmp1 + (ptop*rga)* &
         823     5880600 :                     (h2ostr*zenfac(i)*delta)
         824     5880600 :          uh2o(i,0) = h2ommr(i,1)*pthh2o
         825     5880600 :          uco2(i,0) = zenfac(i)*pthco2
         826     5880600 :          uo2 (i,0) = zenfac(i)*ptho2
         827     6272383 :          uo3 (i,0) = ptho3
         828             : !
         829             : ! End  do i=1,Nday
         830             : !
         831             :    end do
         832             :
         833    10578141 :    do k=1,pver
         834             :
         835   163473741 :       do i=1,Nday
         836             :
         837   152895600 :          pdel      = pflx(i,k+1) - pflx(i,k)
         838   152895600 :          path      = pdel / gravit
         839   152895600 :          ptho2     = o2mmr(i) * path
         840   152895600 :          ptho3     = o3mmr(i,k) * path
         841   152895600 :          pthco2    = sqrco2(i) * path
         842   152895600 :          h2ostr    = sqrt(1.0_r8/h2ommr(i,k))
         843   152895600 :          pthh2o    = (pflx(i,k+1)**2 - pflx(i,k)**2)*tmp1 + pdel*h2ostr*zenfac(i)*tmp2
         844   152895600 :          uh2o(i,k) = h2ommr(i,k)*pthh2o
         845   152895600 :          uco2(i,k) = zenfac(i)*pthco2
         846   152895600 :          uo2 (i,k) = zenfac(i)*ptho2
         847   163081958 :          uo3 (i,k) = ptho3
         848             :
         849             : !
         850             : ! End  do i=1,Nday
         851             : !
         852             :       end do
         853             : !
         854             : ! End  k=1,pver
         855             : !
         856             :    end do
         857             : !
         858             : ! Compute column absorber amounts for the clear sky computation:
         859             : !
         860     6272383 :    do i=1,Nday
         861             :
         862     5880600 :       uth2o(i) = 0.0_r8
         863     5880600 :       uto3(i)  = 0.0_r8
         864     5880600 :       utco2(i) = 0.0_r8
         865     5880600 :       uto2(i)  = 0.0_r8
         866             :
         867   159167983 :       do k=1,pver
         868   152895600 :          uth2o(i) = uth2o(i) + uh2o(i,k)
         869   152895600 :          uto3(i)  = uto3(i)  + uo3(i,k)
         870   152895600 :          utco2(i) = utco2(i) + uco2(i,k)
         871   158776200 :          uto2(i)  = uto2(i)  + uo2(i,k)
         872             : !
         873             : ! End  k=1,pver
         874             : !
         875             :       end do
         876             : !
         877             : ! End  do i=1,Nday
         878             : !
         879             :    end do
         880             : !
         881             : ! Set cloud properties for top (0) layer; so long as tauxcl is zero,
         882             : ! there is no cloud above top of model; the other cloud properties
         883             : ! are arbitrary:
         884             : !
         885     6272383 :       do i=1,Nday
         886             :
         887     5880600 :          tauxcl(i,0)  = 0._r8
         888     5880600 :          wcl(i,0)     = 0.999999_r8
         889     5880600 :          gcl(i,0)     = 0.85_r8
         890     5880600 :          fcl(i,0)     = 0.725_r8
         891     5880600 :          tauxci(i,0)  = 0._r8
         892     5880600 :          wci(i,0)     = 0.999999_r8
         893     5880600 :          gci(i,0)     = 0.85_r8
         894     6272383 :          fci(i,0)     = 0.725_r8
         895             : !
         896             : ! End  do i=1,Nday
         897             : !
         898             :       end do
         899             : !
         900             : ! Begin spectral loop
         901             : !
         902     7835660 :    do ns=1,nswbands
         903             : !
         904             : ! Set cloud extinction optical depth, single scatter albedo,
         905             : ! asymmetry parameter, and forward scattered fraction:
         906             : !
         907     7443877 :       abarli = abarl(indxsl(ns))
         908     7443877 :       bbarli = bbarl(indxsl(ns))
         909     7443877 :       cbarli = cbarl(indxsl(ns))
         910     7443877 :       dbarli = dbarl(indxsl(ns))
         911     7443877 :       ebarli = ebarl(indxsl(ns))
         912     7443877 :       fbarli = fbarl(indxsl(ns))
         913             : !
         914     7443877 :       abarii = abari(indxsl(ns))
         915     7443877 :       bbarii = bbari(indxsl(ns))
         916     7443877 :       cbarii = cbari(indxsl(ns))
         917     7443877 :       dbarii = dbari(indxsl(ns))
         918     7443877 :       ebarii = ebari(indxsl(ns))
         919     7443877 :       fbarii = fbari(indxsl(ns))
         920             : !
         921             : ! adjustfraction within spectral interval to allow for the possibility of
         922             : ! sub-divisions within a particular interval:
         923             : !
         924     7443877 :       psf(ns) = 1.0_r8
         925     7443877 :       if(ph2o(ns)/=0._r8) psf(ns) = psf(ns)*ph2o(ns)
         926     7443877 :       if(pco2(ns)/=0._r8) psf(ns) = psf(ns)*pco2(ns)
         927     7443877 :       if(po2 (ns)/=0._r8) psf(ns) = psf(ns)*po2 (ns)
         928             :
         929   200984679 :       do k=1,pver
         930             :
         931  3106001079 :          do i=1,Nday
         932             :
         933             :             ! liquid
         934             :             ! note that optical properties for liquid valid only
         935             :             ! in range of 4.2 > rel > 16 micron (Slingo 89)
         936  2905016400 :             if ( microp_scheme == 'MG' ) then
         937           0 :                tmp2l = 1._r8 - cbarli - dbarli*min(max(4.2_r8,rel(i,k)),16._r8)
         938           0 :                tmp3l = fbarli*min(max(4.2_r8,rel(i,k)),16._r8)
         939             :             else
         940  2905016400 :                tmp2l = 1._r8 - cbarli - dbarli*rel(i,k)
         941  2905016400 :                tmp3l = fbarli*rel(i,k)
         942             :             endif
         943             :
         944             :             ! ice
         945             :             ! note that optical properties for ice valid only
         946             :             ! in range of 13 > rei > 130 micron (Ebert and Curry 92)
         947  2905016400 :             if ( microp_scheme == 'MG' ) then
         948           0 :                tmp2i = 1._r8 - cbarii - dbarii*min(max(13._r8,rei(i,k)),130._r8)
         949           0 :                tmp3i = fbarii*min(max(13._r8,rei(i,k)),130._r8)
         950             :             else
         951  2905016400 :                tmp2i = 1._r8 - cbarii - dbarii*rei(i,k)
         952  2905016400 :                tmp3i = fbarii*rei(i,k)
         953             :             endif
         954             :
         955  2905016400 :             if (cld(i,k) >= cldmin .and. cld(i,k) >= cldeps) then
         956             :
         957             :                ! liquid
         958   735757748 :                if ( microp_scheme == 'MG' ) then
         959           0 :                   tmp1l = abarli + bbarli/min(max(4.2_r8,rel(i,k)),16._r8)
         960             :                else
         961   735757748 :                   tmp1l = abarli + bbarli/rel(i,k)
         962             :                endif
         963             :
         964             :                ! ice
         965   735757748 :                if ( microp_scheme == 'MG' ) then
         966           0 :                   tmp1i = abarii + bbarii/max(13._r8,min(rei(i,k),130._r8))
         967             :                else
         968   735757748 :                   tmp1i = abarii + bbarii/rei(i,k)
         969             :                endif
         970             :
         971   735757748 :                tauxcl(i,k) = cliqwp(i,k)*tmp1l
         972   735757748 :                tauxci(i,k) = cicewp(i,k)*tmp1i
         973             :
         974   735757748 :                if (ns .eq. idx_sw_diag) then
         975    38724092 :                   tauxcl_out(i,k) = cliqwp(i,k)*tmp1l
         976    38724092 :                   tauxci_out(i,k) = cicewp(i,k)*tmp1i
         977             :                endif
         978             :             else
         979  2169258652 :                tauxcl(i,k) = 0.0_r8
         980  2169258652 :                tauxci(i,k) = 0.0_r8
         981             :             endif
         982             :
         983             :             ! Do not let single scatter albedo be 1.  Delta-eddington solution
         984             :             ! for non-conservative case has different analytic form from solution
         985             :             ! for conservative case, and raddedmx is written for non-conservative case.
         986  2905016400 :             wcl(i,k) = min(tmp2l,.999999_r8)
         987  2905016400 :             gcl(i,k) = ebarli + tmp3l
         988  2905016400 :             fcl(i,k) = gcl(i,k)*gcl(i,k)
         989             :
         990  2905016400 :             wci(i,k) = min(tmp2i,.999999_r8)
         991  2905016400 :             gci(i,k) = ebarii + tmp3i
         992  3098557202 :             fci(i,k) = gci(i,k)*gci(i,k)
         993             :
         994             :          end do ! End do i=1,Nday
         995             :       end do    ! End do k=1,pver
         996             :
         997             :
         998             : !
         999             : ! Set reflectivities for surface based on mid-point wavelength
        1000             : !
        1001     7443877 :       call get_sw_spectral_boundaries(wavmin, wavmax, 'micrometer')
        1002     7443877 :       wavmid(ns) = 0.5_r8*(wavmin(ns) + wavmax(ns))
        1003             : !
        1004             : ! Wavelength less  than 0.7 micro-meter
        1005             : !
        1006     7443877 :       if (wavmid(ns) < 0.7_r8 ) then
        1007    56451447 :          do i=1,Nday
        1008    52925400 :                albdir(i,ns) = asdir(i)
        1009    56451447 :                albdif(i,ns) = asdif(i)
        1010             :          end do
        1011             : !
        1012             : ! Wavelength greater than 0.7 micro-meter
        1013             : !
        1014             :       else
        1015    62723830 :          do i=1,Nday
        1016    58806000 :                albdir(i,ns) = aldir(i)
        1017    62723830 :                albdif(i,ns) = aldif(i)
        1018             :          end do
        1019             :       end if
        1020     7443877 :       trayoslp = raytau(ns)/sslp
        1021             : !
        1022             : ! Layer input properties now completely specified; compute the
        1023             : ! delta-Eddington solution reflectivities and transmissivities
        1024             : ! for each layer
        1025             : !
        1026             :       call raddedmx(coszrs   ,Nday    , &
        1027             :               abh2o(ns),abo3(ns) ,abco2(ns),abo2(ns) , &
        1028             :               uh2o     ,uo3      ,uco2     ,uo2      , &
        1029             :               trayoslp ,pflx     ,ns       , &
        1030             :               tauxcl   ,wcl      ,gcl      ,fcl      , &
        1031             :               tauxci   ,wci      ,gci      ,fci      , &
        1032             :               aer_tau(:,:,ns) ,aer_tau_w(:,:,ns)   ,aer_tau_w_g(:,:,ns) ,aer_tau_w_f(:,:,ns) , &
        1033             :               rdir     ,rdif     ,tdir     ,tdif     ,explay  , &
        1034     7835660 :               rdirc    ,rdifc    ,tdirc    ,tdifc    ,explayc )
        1035             : !
        1036             : ! End spectral loop
        1037             : !
        1038             :    end do
        1039             : !
        1040             : !----------------------------------------------------------------------
        1041             : !
        1042             : ! Solution for max/random cloud overlap.
        1043             : !
        1044             : ! Steps:
        1045             : ! (1. delta-Eddington solution for each layer (called above)
        1046             : !
        1047             : ! (2. The adding method is used to
        1048             : ! compute the reflectivity and transmissivity to direct and diffuse
        1049             : ! radiation from the top and bottom of the atmosphere for each
        1050             : ! cloud configuration.  This calculation is based upon the
        1051             : ! max-random overlap assumption.
        1052             : !
        1053             : ! (3. to solve for the fluxes, combine the
        1054             : ! bulk properties of the atmosphere above/below the region.
        1055             : !
        1056             : ! Index calculations for steps 2-3 are performed outside spectral
        1057             : ! loop to avoid redundant calculations.  Index calculations (with
        1058             : ! application of areamin & nconfgmax conditions) are performed
        1059             : ! first to identify the minimum subset of terms for the configurations
        1060             : ! satisfying the areamin & nconfgmax conditions. This minimum set is
        1061             : ! used to identify the corresponding minimum subset of terms in
        1062             : ! steps 2 and 3.
        1063             : !
        1064     6268528 :    do iconfig = 1, nconfgmax
        1065    94085745 :       ccon(iconfig,0,1:Nday)      = 0
        1066    94085745 :       ccon(iconfig,pverp,1:Nday)  = 0
        1067             :
        1068    94085745 :       icond(iconfig,0,1:Nday)     = iconfig
        1069    94477528 :       iconu(iconfig,pverp,1:Nday) = iconfig
        1070             :    end do
        1071             : !
        1072             : ! Construction of nuniqu/d, istrtu/d, iconu/d using binary tree
        1073             : !
        1074     6272383 :          nuniqd(0,1:Nday) = 1
        1075     6272383 :          nuniqu(pverp,1:Nday) = 1
        1076             :
        1077     6272383 :          istrtd(1,0,1:Nday) = 1
        1078     6272383 :          istrtu(1,pverp,1:Nday) = 1
        1079             :
        1080             :
        1081             : !CSD$ PARALLEL DO PRIVATE( npasses, kx2, mrgn, region_found, k1, k2, kx1, nxs, ksort, asort ) &
        1082             : !CSD$ PRIVATE ( ktmp, atmp, cstr, mstr, nstr, cld0, wstr, nrgn, nconfigm, istr, new_term, xwgt ) &
        1083             : !CSD$ PRIVATE ( j, ptrc, wgtv, km1, nuniq, is0, is1, n0, n1, ptr0, ptr1, kp1, i, irgn ) &
        1084             : !CSD$ PRIVATE ( k, l, iconfig, l0, isn )
        1085     6272383 :    do i=1,Nday
        1086             :
        1087             : !----------------------------------------------------------------------
        1088             : ! INDEX CALCULATIONS FOR MAX OVERLAP
        1089             : !
        1090             : ! The column is divided into sets of adjacent layers, called regions,
        1091             : ! in which the clouds are maximally overlapped.  The clouds are
        1092             : ! randomly overlapped between different regions.  The number of
        1093             : ! regions in a column is set by nmxrgn, and the range of pressures
        1094             : ! included in each region is set by pmxrgn.
        1095             : !
        1096             : ! The following calculations determine the number of unique cloud
        1097             : ! configurations (assuming maximum overlap), called "streams",
        1098             : ! within each region. Each stream consists of a vector of binary
        1099             : ! clouds (either 0 or 100% cloud cover).  Over the depth of the region,
        1100             : ! each stream requires a separate calculation of radiative properties. These
        1101             : ! properties are generated using the adding method from
        1102             : ! the radiative properties for each layer calculated by raddedmx.
        1103             : !
        1104             : ! The upward and downward-propagating streams are treated
        1105             : ! separately.
        1106             : !
        1107             : ! We will refer to a particular configuration of binary clouds
        1108             : ! within a single max-overlapped region as a "stream".  We will
        1109             : ! refer to a particular arrangement of binary clouds over the entire column
        1110             : ! as a "configuration".
        1111             : !
        1112             : ! This section of the code generates the following information:
        1113             : ! (1. nrgn    : the true number of max-overlap regions (need not = nmxrgn)
        1114             : ! (2. nstr    : the number of streams in a region (>=1)
        1115             : ! (3. cstr    : flags for presence of clouds at each layer in each stream
        1116             : ! (4. wstr    : the fractional horizontal area of a grid box covered
        1117             : ! by each stream
        1118             : ! (5. kx1,2   : level indices for top/bottom of each region
        1119             : !
        1120             : ! The max-overlap calculation proceeds in 3 stages:
        1121             : ! (1. compute layer radiative properties in raddedmx.
        1122             : ! (2. combine these properties between layers
        1123             : ! (3. combine properties to compute fluxes at each interface.
        1124             : !
        1125             : ! Most of the indexing information calculated here is used in steps 2-3
        1126             : ! after the call to raddedmx.
        1127             : !
        1128             : ! Initialize indices for layers to be max-overlapped
        1129             : !
        1130             : ! Loop to handle fix in totwgt=0. For original overlap config
        1131             : ! from npasses = 0.
        1132             : !
        1133             :          npasses = 0
        1134           0 :          do
        1135    20810214 :             do irgn = 0, nmxrgn(i)
        1136    20810214 :                kx2(irgn) = 0
        1137             :             end do
        1138             :             mrgn = 0
        1139             : !
        1140             : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
        1141             : !
        1142    14929614 :             do irgn = 1, nmxrgn(i)
        1143             : !
        1144             : ! Calculate min/max layer indices inside region.
        1145             : !
        1146     9049014 :                region_found = .false.
        1147     9049014 :                if (kx2(irgn-1) < pver) then
        1148     9049014 :                   k1 = kx2(irgn-1)+1
        1149     9049014 :                   kx1(irgn) = k1
        1150     9049014 :                   kx2(irgn) = k1-1
        1151    44578071 :                   do k2 = pver, k1, -1
        1152    44578071 :                      if (pmid(i,k2) <= pmxrgn(i,irgn)) then
        1153     9049014 :                         kx2(irgn) = k2
        1154     9049014 :                         mrgn = mrgn+1
        1155             :                         region_found = .true.
        1156             :                         exit
        1157             :                      end if
        1158             :                   end do
        1159             :                else
        1160             :                   exit
        1161             :                endif
        1162             :
        1163     5880600 :                if (region_found) then
        1164             : !
        1165             : ! Sort cloud areas and corresponding level indices.
        1166             : !
        1167     9049014 :                   nxs = 0
        1168             :                   if (cldeps > 0) then
        1169             :                      do k = k1,k2
        1170             :                         if (cld(i,k) >= cldmin .and. cld(i,k) >= cldeps) then
        1171             :                            nxs = nxs+1
        1172             :                            ksort(nxs) = k
        1173             : !
        1174             : ! We need indices for clouds in order of largest to smallest, so
        1175             : ! sort 1-cld in ascending order
        1176             : !
        1177             :                            asort(nxs) = 1.0_r8-(floor(cld(i,k)/cldeps)*cldeps)
        1178             :                         end if
        1179             :                      end do
        1180             :                   else
        1181   161944614 :                      do k = k1,k2
        1182   161944614 :                         if (cld(i,k) >= cldmin) then
        1183    38724092 :                            nxs = nxs+1
        1184    38724092 :                            ksort(nxs) = k
        1185             : !
        1186             : ! We need indices for clouds in order of largest to smallest, so
        1187             : ! sort 1-cld in ascending order
        1188             : !
        1189    38724092 :                            asort(nxs) = 1.0_r8-cld(i,k)
        1190             :                         end if
        1191             :                      end do
        1192             :                   endif
        1193             : !
        1194             : ! If nxs eq 1, no need to sort.
        1195             : ! If nxs eq 2, sort by swapping if necessary
        1196             : ! If nxs ge 3, sort using local sort routine
        1197             : !
        1198     9049014 :                   if (nxs == 2) then
        1199     2052291 :                      if (asort(2) < asort(1)) then
        1200     1550125 :                         ktmp = ksort(1)
        1201     1550125 :                         ksort(1) = ksort(2)
        1202     1550125 :                         ksort(2) = ktmp
        1203             :
        1204     1550125 :                         atmp = asort(1)
        1205     1550125 :                         asort(1) = asort(2)
        1206     1550125 :                         asort(2) = atmp
        1207             :                      endif
        1208     6996723 :                   else if (nxs >= 3) then
        1209     5534586 :                      call quick_sort(asort(1:nxs),ksort(1:nxs))
        1210             :                   endif
        1211             : !
        1212             : ! Construct wstr, cstr, nstr for this region
        1213             : !
        1214   953873116 :                   cstr(k1:k2,1:nxs+1) = 0
        1215             :                   mstr = 1
        1216             :                   cld0 = 0.0_r8
        1217    47773106 :                   do l = 1, nxs
        1218    38724092 :                      if (asort(l) /= cld0) then
        1219    37598493 :                         wstr(mstr,mrgn) = asort(l) - cld0
        1220    37598493 :                         cld0 = asort(l)
        1221    37598493 :                         mstr = mstr + 1
        1222             :                      endif
        1223   207609127 :                      cstr(ksort(l),mstr:nxs+1) = 1
        1224             :                   end do
        1225     9049014 :                   nstr(mrgn) = mstr
        1226     9049014 :                   wstr(mstr,mrgn) = 1.0_r8 - cld0
        1227             : !
        1228             : ! End test of region_found = true
        1229             : !
        1230             :                endif
        1231             : !
        1232             : ! End loop over regions irgn for max-overlap
        1233             : !
        1234             :             end do
        1235     5880600 :             nrgn = mrgn
        1236             : !
        1237             : ! Finish construction of cstr for additional top layer
        1238             : !
        1239    37613302 :             cstr(0,1:nstr(1)) = 0
        1240             : !
        1241             : ! INDEX COMPUTATIONS FOR STEP 2-3
        1242             : ! This section of the code generates the following information:
        1243             : ! (1. totwgt     step 3     total frac. area of configurations satisfying
        1244             : ! areamin & nconfgmax criteria
        1245             : ! (2. wgtv       step 3     frac. area of configurations
        1246             : ! (3. ccon       step 2     binary flag for clouds in each configuration
        1247             : ! (4. nconfig    steps 2-3  number of configurations
        1248             : ! (5. nuniqu/d   step 2     Number of unique cloud configurations for
        1249             : ! up/downwelling rad. between surface/TOA
        1250             : ! and level k
        1251             : ! (6. istrtu/d   step 2     Indices into iconu/d
        1252             : ! (7. iconu/d    step 2     Cloud configurations which are identical
        1253             : ! for up/downwelling rad. between surface/TOA
        1254             : ! and level k
        1255             : !
        1256             : ! Number of configurations (all permutations of streams in each region)
        1257             : !
        1258    14929614 :             nconfigm = product(nstr(1: nrgn))
        1259             : !
        1260             : ! Construction of totwgt, wgtv, ccon, nconfig
        1261             : !
        1262    14929614 :             istr(1: nrgn) = 1
        1263     5880600 :             nconfig(i) = 0
        1264     5880600 :             totwgt(i) = 0.0_r8
        1265     5880600 :             new_term = .true.
        1266    96641319 :             do iconfig = 1, nconfigm
        1267             :                xwgt = 1.0_r8
        1268   281580814 :                do mrgn = 1,  nrgn
        1269   281580814 :                   xwgt = xwgt * wstr(istr(mrgn),mrgn)
        1270             :                end do
        1271    90760719 :                if (xwgt >= areamin) then
        1272    38641753 :                   nconfig(i) = nconfig(i) + 1
        1273    38641753 :                   if (nconfig(i) <= nconfgmax) then
        1274    37842590 :                      j = nconfig(i)
        1275    37842590 :                      ptrc(nconfig(i)) = nconfig(i)
        1276             :                   else
        1277      799163 :                      nconfig(i) = nconfgmax
        1278      799163 :                      if (new_term) then
        1279      666497 :                         min_idx = minloc(wgtv)
        1280      666497 :                         j = min_idx(1)
        1281             :                      endif
        1282      799163 :                      if (wgtv(j) < xwgt) then
        1283      593723 :                         totwgt(i) = totwgt(i) - wgtv(j)
        1284             :                         new_term = .true.
        1285             :                      else
        1286             :                         new_term = .false.
        1287             :                      endif
        1288             :                   endif
        1289    37842590 :                   if (new_term) then
        1290    38436313 :                      wgtv(j) = xwgt
        1291    38436313 :                      totwgt(i) = totwgt(i) + xwgt
        1292   104373107 :                      do mrgn = 1, nrgn
        1293  1103717245 :                         ccon(j,kx1(mrgn):kx2(mrgn),i) = cstr(kx1(mrgn):kx2(mrgn),istr(mrgn))
        1294             :                      end do
        1295             :                   endif
        1296             :                endif
        1297             :
        1298    90760719 :                mrgn =  nrgn
        1299    90760719 :                istr(mrgn) = istr(mrgn) + 1
        1300   113795766 :                do while (istr(mrgn) > nstr(mrgn) .and. mrgn > 1)
        1301    17154447 :                   istr(mrgn) = 1
        1302    17154447 :                   mrgn = mrgn - 1
        1303    17154447 :                   istr(mrgn) = istr(mrgn) + 1
        1304             :                end do
        1305             : !
        1306             : ! End do iconfig = 1, nconfigm
        1307             : !
        1308             :             end do
        1309             : !
        1310             : ! If totwgt(i) = 0 implement maximum overlap and make another pass
        1311             : ! if totwgt(i) = 0 on this second pass then terminate.
        1312             : !
        1313     5880600 :             if (totwgt(i) > 0._r8) then
        1314             :                exit
        1315             :             else
        1316           0 :                npasses = npasses + 1
        1317           0 :                if (npasses >= 2 ) then
        1318           0 :                   write(iulog,*)'RADCSWMX: Maximum overlap of column ','failed'
        1319           0 :                   call endrun('RADCSWMX')
        1320             :                endif
        1321           0 :                nmxrgn(i)=1
        1322           0 :                pmxrgn(i,1)=1.0e30_r8
        1323             :             end if
        1324             : !
        1325             : ! End npasses = 0, do
        1326             : !
        1327             :          end do
        1328             : !
        1329             : ! Finish construction of ccon
        1330             : !
        1331             :
        1332     5880600 :          istrtd(2,0,i) = nconfig(i)+1
        1333     5880600 :          istrtu(2,pverp,i) = nconfig(i)+1
        1334             :
        1335   164656800 :          do k = 1, pverp
        1336   158776200 :             km1 = k-1
        1337   158776200 :             nuniq = 0
        1338   158776200 :             istrtd(1,k,i) = 1
        1339   541164606 :             do l0 = 1, nuniqd(km1,i)
        1340   382388406 :                is0 = istrtd(l0,km1,i)
        1341   382388406 :                is1 = istrtd(l0+1,km1,i)-1
        1342   382388406 :                n0 = 0
        1343   382388406 :                n1 = 0
        1344  1404138336 :                do isn = is0, is1
        1345  1021749930 :                   j = icond(isn,km1,i)
        1346  1404138336 :                   if (ccon(j,k,i) == 0) then
        1347   887016785 :                      n0 = n0 + 1
        1348   887016785 :                      ptr0(n0) = j
        1349             :                   else       ! if (ccon(j,k,i) == 1) then
        1350   134733145 :                      n1 = n1 + 1
        1351   134733145 :                      ptr1(n1) = j
        1352             :                   endif
        1353             :                end do
        1354   382388406 :                if (n0 > 0) then
        1355   328633987 :                   nuniq = nuniq + 1
        1356   328633987 :                   istrtd(nuniq+1,k,i) = istrtd(nuniq,k,i)+n0
        1357  1215650772 :                   icond(istrtd(nuniq,k,i):istrtd(nuniq+1,k,i)-1,k,i) =  ptr0(1:n0)
        1358             :                endif
        1359   541164606 :                if (n1 > 0) then
        1360    85716409 :                   nuniq = nuniq + 1
        1361    85716409 :                   istrtd(nuniq+1,k,i) = istrtd(nuniq,k,i)+n1
        1362   220449554 :                   icond(istrtd(nuniq,k,i):istrtd(nuniq+1,k,i)-1,k,i) =  ptr1(1:n1)
        1363             :                endif
        1364             :             end do
        1365   164656800 :             nuniqd(k,i) = nuniq
        1366             :          end do
        1367             : !
        1368             : !  Find 'transition point' in downward configurations where the number
        1369             : !  of 'configurations' changes from 1.  This is used to optimize the
        1370             : !  construction of the upward configurations.
        1371             : !  Note: k1 == transition point
        1372             : !
        1373             :
        1374    69526800 :          do k = pverp,0,-1
        1375    69526800 :            if ( nuniqd(k,i) == 1) then
        1376             :               k1 = k
        1377             :               exit
        1378             :            end if
        1379             :          end do
        1380             :
        1381    63880230 :          do k = pver, k1+1,-1
        1382    57999630 :             kp1 = k+1
        1383    57999630 :             nuniq = 0
        1384    57999630 :             istrtu(1,k,i) = 1
        1385   253767407 :             do l0 = 1, nuniqu(kp1,i)
        1386   195767777 :                is0 = istrtu(l0,kp1,i)
        1387   195767777 :                is1 = istrtu(l0+1,kp1,i)-1
        1388   195767777 :                n0 = 0
        1389   195767777 :                n1 = 0
        1390   643879007 :                do isn = is0, is1
        1391   448111230 :                   j = iconu(isn,kp1,i)
        1392   643879007 :                   if (ccon(j,k,i) == 0) then
        1393   314719835 :                      n0 = n0 + 1
        1394   314719835 :                      ptr0(n0) = j
        1395             :                   else       ! if (ccon(j,k,i) == 1) then
        1396   133391395 :                      n1 = n1 + 1
        1397   133391395 :                      ptr1(n1) = j
        1398             :                   endif
        1399             :                end do
        1400   195767777 :                if (n0 > 0) then
        1401   145872508 :                   nuniq = nuniq + 1
        1402   145872508 :                   istrtu(nuniq+1,k,i) = istrtu(nuniq,k,i)+n0
        1403   460592343 :                   iconu(istrtu(nuniq,k,i):istrtu(nuniq+1,k,i)-1,k,i) =  ptr0(1:n0)
        1404             :                endif
        1405   253767407 :                if (n1 > 0) then
        1406    81857259 :                   nuniq = nuniq + 1
        1407    81857259 :                   istrtu(nuniq+1,k,i) = istrtu(nuniq,k,i)+n1
        1408   215248654 :                   iconu(istrtu(nuniq,k,i):istrtu(nuniq+1,k,i)-1,k,i) = ptr1(1:n1)
        1409             :                endif
        1410             :             end do
        1411    63880230 :             nuniqu(k,i) = nuniq
        1412             :          end do
        1413             : !
        1414             : !  Copy identical configurations from 'transition point' to surface.
        1415             : !
        1416     5880600 :          k1 = min(pverp-1,k1)
        1417     5880600 :          nuniq = nuniqu(k1+1,i)
        1418   106657170 :          do k = k1,0,-1
        1419   100776570 :             nuniqu(k,i) = nuniq
        1420   674415270 :             iconu(1:nuniq,k,i) = iconu(1:nuniq,k1+1,i)
        1421   781072440 :             istrtu(1:nuniq+1,k,i) = istrtu(1:nuniq+1,k1+1,i)
        1422             :          end do
        1423             :
        1424    44114973 :          v_wgtv(1:nconfig(i),i) = wgtv(1:nconfig(i))
        1425             :
        1426             : !
        1427             : !----------------------------------------------------------------------
        1428             : ! End of index calculations
        1429             : !----------------------------------------------------------------------
        1430             : !
        1431             : ! End do i=1,Nday
        1432             : !
        1433             :    end do
        1434             : !CSD$ END PARALLEL
        1435             :
        1436             : !----------------------------------------------------------------------
        1437             : ! Start of flux calculations
        1438             : !----------------------------------------------------------------------
        1439             : !
        1440             : ! Initialize spectrally integrated totals:
        1441             : !
        1442   169746124 :          totfld (1:Nday,0:pver) = 0.0_r8
        1443   169746124 :          totfldc(1:Nday,0:pver) = 0.0_r8
        1444   169746124 :          fswup  (1:Nday,0:pver) = 0.0_r8
        1445   169746124 :          fswdn  (1:Nday,0:pver) = 0.0_r8
        1446             :
        1447     6272383 :          fswup (1:Nday,pverp)  = 0.0_r8
        1448     6272383 :          fswdn (1:Nday,pverp)  = 0.0_r8
        1449             : !
        1450             : ! Start spectral interval
        1451             : !
        1452             : !old   do ns = 1,nswbands
        1453             : !old     wgtint = nirwgt(ns)
        1454             :
        1455     6272383 :      do i=1,Nday
        1456             :
        1457             : !----------------------------------------------------------------------
        1458             : ! STEP 2
        1459             : !
        1460             : !
        1461             : ! Apply adding method to solve for radiative properties
        1462             : !
        1463             : ! first initialize the bulk properties at toa
        1464             : !
        1465             :
        1466             : ! nswbands, 0:pverp, nconfgmax, pcols
        1467             :
        1468   762732400 :             rdndif(:,0,1:nconfig(i),i) = 0.0_r8
        1469   762732400 :             exptdn(:,0,1:nconfig(i),i) = 1.0_r8
        1470   763124183 :             tdntot(:,0,1:nconfig(i),i) = 1.0_r8
        1471             : !
        1472             : ! End do i=1,Nday
        1473             : !
        1474             :      end do
        1475             : !
        1476             : ! solve for properties involving downward propagation of radiation.
        1477             : ! the bulk properties are:
        1478             : !
        1479             : ! (1. exptdn   sol. beam dwn. trans from layers above
        1480             : ! (2. rdndif   ref to dif rad for layers above
        1481             : ! (3. tdntot   total trans for layers above
        1482             : !
        1483             :
        1484             : !CSD$ PARALLEL DO PRIVATE( km1, is0, is1, j, jj, Ttdif, Trdif, Trdir, Ttdir, Texplay ) &
        1485             : !CSD$ PRIVATE( xexpt, xrdnd, tdnmexp,  ytdnd, yrdnd, rdenom, rdirexp, zexpt, zrdnd, ztdnt ) &
        1486             : !CSD$ PRIVATE( i, k, l0, ns, isn )
        1487     6272383 :          do i = 1, Nday
        1488   165048583 :             do k = 1, pverp
        1489   158776200 :                km1 = k - 1
        1490   547045206 :                do l0 = 1, nuniqd(km1,i)
        1491   382388406 :                   is0 = istrtd(l0,km1,i)
        1492   382388406 :                   is1 = istrtd(l0+1,km1,i)-1
        1493             :
        1494   382388406 :                   j = icond(is0,km1,i)
        1495             :
        1496             : !
        1497             : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
        1498             : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
        1499             : !
        1500   382388406 :                   if ( ccon(j,km1,i) == 1 ) then
        1501  1714328180 :                      Ttdif(:) = tdif(:,i,km1)
        1502  1714328180 :                      Trdif(:) = rdif(:,i,km1)
        1503  1714328180 :                      Trdir(:) = rdir(:,i,km1)
        1504  1714328180 :                      Ttdir(:) = tdir(:,i,km1)
        1505  1714328180 :                      Texplay(:) = explay(:,i,km1)
        1506             :                   else
        1507  5933439940 :                      Ttdif(:) = tdifc(:,i,km1)
        1508  5933439940 :                      Trdif(:) = rdifc(:,i,km1)
        1509  5933439940 :                      Trdir(:) = rdirc(:,i,km1)
        1510  5933439940 :                      Ttdir(:) = tdirc(:,i,km1)
        1511  5933439940 :                      Texplay(:) = explayc(:,i,km1)
        1512             :                   end if
        1513             :
        1514  7647768120 :                   do ns = 1, nswbands
        1515  7265379714 :                   xexpt   = exptdn(ns,km1,j,i)
        1516  7265379714 :                   xrdnd   = rdndif(ns,km1,j,i)
        1517  7265379714 :                   tdnmexp = tdntot(ns,km1,j,i) - xexpt
        1518             :
        1519  7265379714 :                   ytdnd = Ttdif(ns)
        1520  7265379714 :                   yrdnd = Trdif(ns)
        1521             :
        1522  7265379714 :                   rdenom  = 1._r8/(1._r8-yrdnd*xrdnd)
        1523  7265379714 :                   rdirexp = Trdir(ns)*xexpt
        1524             :
        1525  7265379714 :                   zexpt = xexpt * Texplay(ns)
        1526  7265379714 :                   zrdnd = yrdnd + xrdnd*(ytdnd**2)*rdenom
        1527             :                   ztdnt = xexpt*Ttdir(ns) + ytdnd* &
        1528  7265379714 :                           (tdnmexp + xrdnd*rdirexp)*rdenom
        1529             :
        1530  7265379714 :                   exptdn(ns,k,j,i) = zexpt
        1531  7265379714 :                   rdndif(ns,k,j,i) = zrdnd
        1532  7647768120 :                   tdntot(ns,k,j,i) = ztdnt
        1533             :                   end do ! ns = 1, nswbands
        1534             : !
        1535             : ! If 2 or more configurations share identical properties at a given level k,
        1536             : ! the properties (at level k) are computed once and copied to
        1537             : ! all the configurations for efficiency.
        1538             : !
        1539  1180526130 :                   do isn = is0+1, is1
        1540   639361524 :                      jj = icond(isn,km1,i)
        1541 12787230480 :                      exptdn(:,k,jj,i) = exptdn(:,k,j,i)
        1542 12787230480 :                      rdndif(:,k,jj,i) = rdndif(:,k,j,i)
        1543 13169618886 :                      tdntot(:,k,jj,i) = tdntot(:,k,j,i)
        1544             :                   end do
        1545             :
        1546             : !
        1547             : ! end do l0 = 1, nuniqd(k,i)
        1548             : !
        1549             :                end do
        1550             : !
        1551             : ! end do k = 1, pverp
        1552             : !
        1553             :             end do
        1554             : !
        1555             : ! end do i = 1, Nday
        1556             : !
        1557             :          end do
        1558             : !CSD$ END PARALLEL
        1559             : !
        1560             : ! Solve for properties involving upward propagation of radiation.
        1561             : ! The bulk properties are:
        1562             : !
        1563             : ! (1. rupdif   Ref to dif rad for layers below
        1564             : ! (2. rupdir   Ref to dir rad for layers below
        1565             : !
        1566             : ! Specify surface boundary conditions (surface albedos)
        1567             : !
        1568             :
        1569             : ! nswbands, 0:pverp, nconfgmax, pcols
        1570      391783 :    rupdir = 0._r8
        1571      391783 :    rupdif = 0._r8
        1572     6272383 :    do i = 1, Nday
        1573   118003783 :       do ns = 1, nswbands
        1574   830740610 :          rupdir(ns,pverp,1:nconfig(i),i) = albdir(i,ns)
        1575   836621210 :          rupdif(ns,pverp,1:nconfig(i),i) = albdif(i,ns)
        1576             :       end do
        1577             :    end do
        1578             :
        1579     6272383 :          do i = 1, Nday
        1580   165048583 :             do k = pver, 0, -1
        1581   966025267 :                do l0 = 1, nuniqu(k,i)
        1582   801368467 :                   is0 = istrtu(l0,k,i)
        1583   801368467 :                   is1 = istrtu(l0+1,k,i)-1
        1584             :
        1585   801368467 :                   j = iconu(is0,k,i)
        1586             :
        1587             : !
        1588             : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
        1589             : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
        1590             : !
        1591   801368467 :                   if ( ccon(j,k,i) == 1 ) then
        1592  1663980180 :                      Ttdif(:) = tdif(:,i,k)
        1593  1663980180 :                      Trdif(:) = rdif(:,i,k)
        1594  1663980180 :                      Trdir(:) = rdir(:,i,k)
        1595  1663980180 :                      Ttdir(:) = tdir(:,i,k)
        1596  1663980180 :                      Texplay(:) = explay(:,i,k)
        1597             :                   else
        1598 14363389160 :                      Ttdif(:) = tdifc(:,i,k)
        1599 14363389160 :                      Trdif(:) = rdifc(:,i,k)
        1600 14363389160 :                      Trdir(:) = rdirc(:,i,k)
        1601 14363389160 :                      Ttdir(:) = tdirc(:,i,k)
        1602 14363389160 :                      Texplay(:) = explayc(:,i,k)
        1603             :                   end if
        1604             :
        1605 16027369340 :                   do ns = 1, nswbands
        1606 15226000873 :                   xrupd = rupdif(ns,k+1,j,i)
        1607 15226000873 :                   xrups = rupdir(ns,k+1,j,i)
        1608             :
        1609             : !
        1610             : ! If cloud in layer, use cloudy layer radiative properties (ccon == 1)
        1611             : ! If clear layer, use clear-sky layer radiative properties (ccon /= 1)
        1612             : !
        1613 15226000873 :                   yexpt = Texplay(ns)
        1614 15226000873 :                   yrupd = Trdif(ns)
        1615 15226000873 :                   ytupd = Ttdif(ns)
        1616             :
        1617 15226000873 :                   rdenom  = 1._r8/( 1._r8 - yrupd*xrupd)
        1618 15226000873 :                   tdnmexp = (Ttdir(ns)-yexpt)
        1619 15226000873 :                   rdirexp = xrups*yexpt
        1620             :
        1621 15226000873 :                   zrupd = yrupd + xrupd*(ytupd**2)*rdenom
        1622 15226000873 :                   zrups = Trdir(ns) + ytupd*(rdirexp + xrupd*tdnmexp)*rdenom
        1623             :
        1624 15226000873 :                   rupdif(ns,k,j,i) = zrupd
        1625 16027369340 :                   rupdir(ns,k,j,i) = zrups
        1626             :                   end do ! ns = 1, nswbands
        1627             : !
        1628             : ! If 2 or more configurations share identical properties at a given level k,
        1629             : ! the properties (at level k) are computed once and copied to
        1630             : ! all the configurations for efficiency.
        1631             : !
        1632  1180526130 :                   do isn = is0+1, is1
        1633   220381463 :                      jj = iconu(isn,k,i)
        1634  4407629260 :                      rupdif(:,k,jj,i) = rupdif(:,k,j,i)
        1635  5208997727 :                      rupdir(:,k,jj,i) = rupdir(:,k,j,i)
        1636             :                   end do
        1637             :
        1638             : !
        1639             : ! end do l0 = 1, nuniqu(k,i)
        1640             : !
        1641             :                end do
        1642             : !
        1643             : ! end do k = pver,0,-1
        1644             : !
        1645             :             end do
        1646             : !
        1647             : ! end do i = 1, Nday
        1648             : !
        1649             :          end do
        1650             : !
        1651             : !----------------------------------------------------------------------
        1652             : !
        1653             : ! STEP 3
        1654             : !
        1655             : ! Compute up and down fluxes for each interface k.  This requires
        1656             : ! adding up the contributions from all possible permutations
        1657             : ! of streams in all max-overlap regions, weighted by the
        1658             : ! product of the fractional areas of the streams in each region
        1659             : ! (the random overlap assumption).  The adding principle has been
        1660             : ! used in step 2 to combine the bulk radiative properties
        1661             : ! above and below the interface.
        1662             : !
        1663             :
        1664             : !
        1665             : ! Initialize the fluxes
        1666             : !
        1667      391783 :             fluxup = 0.0_r8
        1668      391783 :             fluxdn = 0.0_r8
        1669             :
        1670     6272383 :             do i = 1, Nday
        1671    44114973 :             do iconfig = 1, nconfig(i)
        1672    37842590 :                xwgt = v_wgtv(iconfig,i)
        1673             :
        1674  1103315710 :                do k = 0, pverp
        1675 21229692990 :                   do ns = 1, nswbands
        1676 20132257880 :                   xexpt = exptdn(ns,k,iconfig,i)
        1677 20132257880 :                   xtdnt = tdntot(ns,k,iconfig,i)
        1678 20132257880 :                   xrdnd = rdndif(ns,k,iconfig,i)
        1679 20132257880 :                   xrupd = rupdif(ns,k,iconfig,i)
        1680 20132257880 :                   xrups = rupdir(ns,k,iconfig,i)
        1681             : !
        1682             : ! Flux computation
        1683             : !
        1684 20132257880 :                   rdenom = 1._r8/(1._r8 - xrdnd * xrupd)
        1685             :
        1686             :                   fluxup(ns,k,i) = fluxup(ns,k,i) + xwgt *  &
        1687 20132257880 :                               ((xexpt * xrups + (xtdnt - xexpt) * xrupd) * rdenom)
        1688             :                   fluxdn(ns,k,i) = fluxdn(ns,k,i) + xwgt *  &
        1689 21191850400 :                               (xexpt + (xtdnt - xexpt + xexpt * xrups * xrdnd) * rdenom)
        1690             :                   end do ! do ns = 1, nswbands
        1691             :                end do
        1692             : !
        1693             : ! End do iconfig = 1, nconfig(i)
        1694             : !
        1695             :             end do
        1696             : !
        1697             : ! End do iconfig = 1, Nday
        1698             : !
        1699             :             end do
        1700             : !
        1701             : ! Normalize by total area covered by cloud configurations included
        1702             : ! in solution
        1703             : !
        1704     6272383 :             do i = 1, Nday
        1705   170929183 :             do k = 0, pverp
        1706  3299016600 :             do ns = 1, nswbands
        1707  3128479200 :                fluxup(ns,k,i)=fluxup(ns,k,i) / totwgt(i)
        1708  3293136000 :                fluxdn(ns,k,i)=fluxdn(ns,k,i) / totwgt(i)
        1709             :             end do ! do i = 1, nday
        1710             :             end do ! do k = 0, pverp
        1711             :             end do ! do i = 1, nday
        1712             :
        1713             :
        1714             : !
        1715             : ! Initialize the direct-beam flux at surface
        1716             : !
        1717   118003783 :             wexptdn(:,1:Nday) = 0.0_r8
        1718             :
        1719     7835660 :    do ns = 1,nswbands
        1720     7443877 :      wgtint = nirwgt(ns)
        1721             :
        1722             :
        1723   119175277 :             do i=1,Nday
        1724   838184487 :             do iconfig = 1, nconfig(i)
        1725             : !
        1726             : ! Note: exptdn can be directly indexed by iconfig at k=pverp.
        1727             : !
        1728   830740610 :                wexptdn(ns,i) =  wexptdn(ns,i) + v_wgtv(iconfig,i) * exptdn(ns,pverp,iconfig,i)
        1729             :             end do
        1730             :             end do
        1731             :
        1732   119175277 :             do i=1,Nday
        1733   111731400 :                wexptdn(ns,i) = wexptdn(ns,i) / totwgt(i)
        1734             : !
        1735             : ! Monochromatic computation completed; accumulate in totals
        1736             : !
        1737   111731400 :             if ( do_spctrl_scaling ) then
        1738           0 :                solflx(i)   = solin(i)*frcsol(ns)*psf(ns)*sfac(ns)
        1739             :             else
        1740   111731400 :                solflx(i)   = solin(i)*frcsol(ns)*psf(ns)
        1741             :             endif
        1742   111731400 :             fsnt(i)  = fsnt(i) + solflx(i)*(fluxdn(ns,1,i) - fluxup(ns,1,i))
        1743   111731400 :             fsntoa(i)= fsntoa(i) + solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
        1744   111731400 :             fsutoa(i)= fsutoa(i) + solflx(i)*(fluxup(ns,0,i))
        1745   111731400 :             fsns(i)  = fsns(i) + solflx(i)*(fluxdn(ns,pverp,i)-fluxup(ns,pverp,i))
        1746   111731400 :             fsdtoa(i)  = fsdtoa(i) + solflx(i)
        1747   111731400 :             fswup(i,0) = fswup(i,0) + solflx(i)*fluxup(ns,0,i)
        1748   111731400 :             fswdn(i,0) = fswdn(i,0) + solflx(i)*fluxdn(ns,0,i)
        1749             : !
        1750             : ! Down spectral fluxes need to be in mks; thus the .001 conversion factors
        1751             : !
        1752   111731400 :             if (wavmid(ns) < 0.7_r8) then
        1753    52925400 :                sols(i)  = sols(i) + wexptdn(ns,i)*solflx(i)*0.001_r8
        1754    52925400 :                solsd(i) = solsd(i)+(fluxdn(ns,pverp,i)-wexptdn(ns,i))*solflx(i)*0.001_r8
        1755             :             else
        1756    58806000 :                soll(i)  = soll(i) + wexptdn(ns,i)*solflx(i)*0.001_r8
        1757    58806000 :                solld(i) = solld(i)+(fluxdn(ns,pverp,i)-wexptdn(ns,i))*solflx(i)*0.001_r8
        1758    58806000 :                fsnrtoaq(i) = fsnrtoaq(i) + solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
        1759             :             end if
        1760   119175277 :             fsnirtoa(i) = fsnirtoa(i) + wgtint*solflx(i)*(fluxdn(ns,0,i) - fluxup(ns,0,i))
        1761             :
        1762             : !
        1763             : ! End do i=1,Nday
        1764             : !
        1765             :    end do
        1766             :
        1767             :
        1768   208428556 :             do k=0,pver
        1769  3225176356 :             do i=1,Nday
        1770             : !
        1771             : ! Compute flux divergence in each layer using the interface up and down
        1772             : ! fluxes:
        1773             : !
        1774  3016747800 :                kp1 = k+1
        1775  3016747800 :                flxdiv = (fluxdn(ns,k,i) - fluxdn(ns,kp1,i)) + (fluxup(ns,kp1,i) - fluxup(ns,k,i))
        1776  3016747800 :                totfld(i,k)  = totfld(i,k)  + solflx(i)*flxdiv
        1777  3016747800 :                fswdn(i,kp1) = fswdn(i,kp1) + solflx(i)*fluxdn(ns,kp1,i)
        1778  3016747800 :                fswup(i,kp1) = fswup(i,kp1) + solflx(i)*fluxup(ns,kp1,i)
        1779  3016747800 :                fns(i,kp1)   = fswdn(i,kp1) - fswup(i,kp1)
        1780  3217732479 :                if (single_column.and.scm_crm_mode) then
        1781           0 :                   fus(i,kp1)=fswup(i,kp1)
        1782           0 :                   fds(i,kp1)=fswdn(i,kp1)
        1783             :                endif
        1784             :             end do
        1785             :             end do
        1786             : !
        1787             : ! Perform clear-sky calculation
        1788             : !
        1789             :
        1790   119175277 :             exptdnc(1:Nday,0) =   1.0_r8
        1791   119175277 :             rdndifc(1:Nday,0) =   0.0_r8
        1792   119175277 :             tdntotc(1:Nday,0) =   1.0_r8
        1793   119175277 :             rupdirc(1:Nday,pverp) = albdir(1:Nday,ns)
        1794   119175277 :             rupdifc(1:Nday,pverp) = albdif(1:Nday,ns)
        1795             :
        1796             :
        1797   208428556 :             do k = 1, pverp
        1798  3225176356 :             do i=1,Nday
        1799  3016747800 :                km1 = k - 1
        1800  3016747800 :                xexpt = exptdnc(i,km1)
        1801  3016747800 :                xrdnd = rdndifc(i,km1)
        1802  3016747800 :                yrdnd = rdifc(ns,i,km1)
        1803  3016747800 :                ytdnd = tdifc(ns,i,km1)
        1804             :
        1805  3016747800 :                exptdnc(i,k) = xexpt*explayc(ns,i,km1)
        1806             :
        1807  3016747800 :                rdenom  = 1._r8/(1._r8 - yrdnd*xrdnd)
        1808  3016747800 :                rdirexp = rdirc(ns,i,km1)*xexpt
        1809  3016747800 :                tdnmexp = tdntotc(i,km1) - xexpt
        1810             :
        1811             :                tdntotc(i,k) = xexpt*tdirc(ns,i,km1) + ytdnd*(tdnmexp + xrdnd*rdirexp)* &
        1812  3016747800 :                                 rdenom
        1813  3217732479 :                rdndifc(i,k) = yrdnd + xrdnd*(ytdnd**2)*rdenom
        1814             : !
        1815             : ! End do i=1,Nday
        1816             : !
        1817             :             end do
        1818             :             end do
        1819             :
        1820   208428556 :             do k=pver,0,-1
        1821  3225176356 :             do i=1,Nday
        1822  3016747800 :                xrupd = rupdifc(i,k+1)
        1823  3016747800 :                yexpt = explayc(ns,i,k)
        1824  3016747800 :                yrupd = rdifc(ns,i,k)
        1825  3016747800 :                ytupd = tdifc(ns,i,k)
        1826             :
        1827  3016747800 :                rdenom = 1._r8/( 1._r8 - yrupd*xrupd)
        1828             :
        1829  3016747800 :                rupdirc(i,k) = rdirc(ns,i,k) + ytupd*(rupdirc(i,k+1)*yexpt + &
        1830  3016747800 :                             xrupd*(tdirc(ns,i,k)-yexpt))*rdenom
        1831  3217732479 :                rupdifc(i,k) = yrupd + xrupd*ytupd**2*rdenom
        1832             : !
        1833             : ! End do i=1,Nday
        1834             : !
        1835             :             end do
        1836             :             end do
        1837             :
        1838   215872433 :             do k=0,pverp
        1839  3344351633 :             do i=1,Nday
        1840  3128479200 :                rdenom    = 1._r8/(1._r8 - rdndifc(i,k)*rupdifc(i,k))
        1841             :                fluxup(ns,k,i) = (exptdnc(i,k)*rupdirc(i,k) + (tdntotc(i,k)-exptdnc(i,k))*rupdifc(i,k))* &
        1842  3128479200 :                            rdenom
        1843             :                fluxdn(ns,k,i) = exptdnc(i,k) + &
        1844             :                            (tdntotc(i,k) - exptdnc(i,k) + exptdnc(i,k)*rupdirc(i,k)*rdndifc(i,k))* &
        1845  3336907756 :                            rdenom
        1846             : !
        1847             : ! End do i=1,Nday
        1848             : !
        1849             :             end do
        1850             :             end do
        1851             :
        1852   208428556 :             do k=0,pver
        1853  3225176356 :             do i=1,Nday
        1854             : !
        1855             : ! Compute flux divergence in each layer using the interface up and down
        1856             : ! fluxes:
        1857             : !
        1858  3016747800 :                kp1 = k+1
        1859  3016747800 :                flxdiv        = (fluxdn(ns,k,i) - fluxdn(ns,kp1,i)) + (fluxup(ns,kp1,i) - fluxup(ns,k,i))
        1860  3217732479 :                totfldc(i,k)  = totfldc(i,k)  + solflx(i)*flxdiv
        1861             :             end do
        1862             :             end do
        1863             :
        1864   119175277 :             do i=1,Nday
        1865   111731400 :             fsntc(i)    = fsntc(i)+solflx(i)*(fluxdn(ns,1,i)-fluxup(ns,1,i))
        1866   111731400 :             fsntoac(i)  = fsntoac(i)+solflx(i)*(fluxdn(ns,0,i)-fluxup(ns,0,i))
        1867   111731400 :             fsnsc(i)    = fsnsc(i)+solflx(i)*(fluxdn(ns,pverp,i)-fluxup(ns,pverp,i))
        1868   111731400 :             fsdsc(i)    = fsdsc(i)+solflx(i)*(fluxdn(ns,pverp,i))
        1869   111731400 :             fsnrtoac(i) = fsnrtoac(i)+wgtint*solflx(i)*(fluxdn(ns,0,i)-fluxup(ns,0,i))
        1870   119175277 :             if (single_column.and.scm_crm_mode) then
        1871           0 :                do k = 1,pverp
        1872           0 :                   fusc(i,k)=fusc(i,k) + solflx(i) * fluxup(ns,k,i)
        1873           0 :                   fdsc(i,k)=fdsc(i,k) + solflx(i) * fluxdn(ns,k,i)
        1874             :                enddo
        1875             :             endif
        1876             : !
        1877             : ! End do i=1,Nday
        1878             : !
        1879             :             end do
        1880             :
        1881   208820339 :             do k = 1,pverp
        1882  3225176356 :             do i=1,Nday
        1883  3217732479 :                fcns(i,k)=fcns(i,k) + solflx(i)*(fluxdn(ns,k,i)-fluxup(ns,k,i))
        1884             :             enddo
        1885             :             enddo
        1886             : !
        1887             : ! End of clear sky calculation
        1888             : !
        1889             :
        1890             : !
        1891             : ! End of spectral interval loop
        1892             : !
        1893             :          end do
        1894             :
        1895     6272383 :    do i=1,Nday
        1896             :
        1897             : !
        1898             : ! Compute solar heating rate (J/kg/s)
        1899             : !
        1900   158776200 :          do k=1,pver
        1901   152895600 :             qrs(i,k) = -1.E-4_r8*gravit*totfld(i,k)/(pint(i,k) - pint(i,k+1))
        1902   158776200 :             qrsc(i,k) = -1.E-4_r8*gravit*totfldc(i,k)/(pint(i,k) - pint(i,k+1))
        1903             :          end do
        1904             : !
        1905             : ! Set the downwelling flux at the surface
        1906             : !
        1907     6272383 :          fsds(i) = fswdn(i,pverp)
        1908             : !
        1909             : ! End do i=1,Nday
        1910             : !
        1911             :    end do
        1912             : !
        1913             : ! Rearrange output arrays.
        1914             : !
        1915             : ! intent(inout)
        1916             : !
        1917      391783 :    call ExpDayNite(pmxrgn,      Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1918      391783 :    call ExpDayNite(nmxrgn,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1919             : !
        1920             : ! intent(out)
        1921             : !
        1922      391783 :    call ExpDayNite(solin,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1923      391783 :    call ExpDayNite(qrs,         Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
        1924      391783 :    call ExpDayNite(qrsc,        Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
        1925      391783 :    call ExpDayNite(fns,         Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1926      391783 :    call ExpDayNite(fcns,        Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1927      391783 :    call ExpDayNite(fsns,        Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1928      391783 :    call ExpDayNite(fsnt,        Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1929      391783 :    call ExpDayNite(fsntoa,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1930      391783 :    call ExpDayNite(fsutoa,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1931      391783 :    call ExpDayNite(fsds,        Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1932      391783 :    call ExpDayNite(fsnsc,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1933      391783 :    call ExpDayNite(fsdsc,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1934      391783 :    call ExpDayNite(fsntc,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1935      391783 :    call ExpDayNite(fsdtoa,      Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1936      391783 :    call ExpDayNite(fsntoac,     Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1937      391783 :    call ExpDayNite(sols,        Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1938      391783 :    call ExpDayNite(soll,        Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1939      391783 :    call ExpDayNite(solsd,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1940      391783 :    call ExpDayNite(solld,       Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1941      391783 :    call ExpDayNite(fsnirtoa,    Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1942      391783 :    call ExpDayNite(fsnrtoac,    Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1943      391783 :    call ExpDayNite(fsnrtoaq,    Nday, IdxDay, Nnite, IdxNite, 1, pcols)
        1944      391783 :    call ExpDayNite(tauxcl_out,  Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
        1945      391783 :    call ExpDayNite(tauxci_out,  Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pver)
        1946             :
        1947             : !  these outfld calls don't work for spmd only outfield in scm mode (nonspmd)
        1948      391783 :    if (single_column.and.scm_crm_mode) then
        1949             :    ! Following outputs added for CRM
        1950           0 :       call ExpDayNite(fus,      Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1951           0 :       call ExpDayNite(fds,      Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1952           0 :       call ExpDayNite(fusc,     Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1953           0 :       call ExpDayNite(fdsc,     Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp)
        1954           0 :       call outfld('FUS     ', fus*1.e-3_r8,  pcols, lchnk)
        1955           0 :       call outfld('FDS     ', fds*1.e-3_r8,  pcols, lchnk)
        1956           0 :       call outfld('FUSC    ', fusc*1.e-3_r8, pcols, lchnk)
        1957           0 :       call outfld('FDSC    ', fdsc*1.e-3_r8, pcols, lchnk)
        1958             :    endif
        1959             : !  write(iulog, '(a, x, i3)') 'radcswmx : exiting, chunk identifier', lchnk
        1960             :
        1961             :    return
        1962      749232 : end subroutine radcswmx
        1963             :
        1964             : !-------------------------------------------------------------------------------
        1965             :
        1966     7443877 : subroutine raddedmx(coszrs  ,ndayc   ,abh2o   , &
        1967             :                     abo3    ,abco2   ,abo2    ,uh2o    ,uo3     , &
        1968             :                     uco2    ,uo2     ,trayoslp,pflx    ,ns      , &
        1969             :                     tauxcl  ,wcl     ,gcl     ,fcl     ,tauxci  , &
        1970             :                     wci     ,gci     ,fci     ,aer_tau ,aer_tau_w, &
        1971             :                     aer_tau_w_g, aer_tau_w_f      ,rdir    ,rdif    ,tdir    , &
        1972             :                     tdif    ,explay  ,rdirc   ,rdifc   ,tdirc   , &
        1973             :                     tdifc   ,explayc )
        1974             : !-----------------------------------------------------------------------
        1975             : !
        1976             : ! Purpose:
        1977             : ! Computes layer reflectivities and transmissivities, from the top down
        1978             : ! to the surface using the delta-Eddington solutions for each layer
        1979             : !
        1980             : ! Method:
        1981             : ! For more details , see Briegleb, Bruce P., 1992: Delta-Eddington
        1982             : ! Approximation for Solar Radiation in the NCAR Community Climate Model,
        1983             : ! Journal of Geophysical Research, Vol 97, D7, pp7603-7612).
        1984             : !
        1985             : ! Modified for maximum/random cloud overlap by Bill Collins and John
        1986             : !    Truesdale
        1987             : !
        1988             : ! Author: Bill Collins
        1989             : !
        1990             : !-----------------------------------------------------------------------
        1991             :
        1992             : !
        1993             : ! Minimum total transmission below which no layer computation are done:
        1994             : !
        1995             :    real(r8) trmin                ! Minimum total transmission allowed
        1996             :    real(r8) wray                 ! Rayleigh single scatter albedo
        1997             :    real(r8) gray                 ! Rayleigh asymmetry parameter
        1998             :    real(r8) fray                 ! Rayleigh forward scattered fraction
        1999             :
        2000             :    parameter (trmin = 1.e-3_r8)
        2001             :    parameter (wray = 0.999999_r8)
        2002             :    parameter (gray = 0.0_r8)
        2003             :    parameter (fray = 0.1_r8)
        2004             : !
        2005             : !------------------------------Arguments--------------------------------
        2006             : !
        2007             : ! Input arguments
        2008             : !
        2009             :    real(r8), intent(in) :: coszrs(pcols)        ! Cosine zenith angle
        2010             :    real(r8), intent(in) :: trayoslp             ! Tray/sslp
        2011             :    real(r8), intent(in) :: pflx(pcols,0:pverp)  ! Interface pressure
        2012             :    real(r8), intent(in) :: abh2o                ! Absorption coefficiant for h2o
        2013             :    real(r8), intent(in) :: abo3                 ! Absorption coefficiant for o3
        2014             :    real(r8), intent(in) :: abco2                ! Absorption coefficiant for co2
        2015             :    real(r8), intent(in) :: abo2                 ! Absorption coefficiant for o2
        2016             :    real(r8), intent(in) :: uh2o(pcols,0:pver)   ! Layer absorber amount of h2o
        2017             :    real(r8), intent(in) :: uo3(pcols,0:pver)    ! Layer absorber amount of  o3
        2018             :    real(r8), intent(in) :: uco2(pcols,0:pver)   ! Layer absorber amount of co2
        2019             :    real(r8), intent(in) :: uo2(pcols,0:pver)    ! Layer absorber amount of  o2
        2020             :    real(r8), intent(in) :: tauxcl(pcols,0:pver) ! Cloud extinction optical depth (liquid)
        2021             :    real(r8), intent(in) :: wcl(pcols,0:pver)    ! Cloud single scattering albedo (liquid)
        2022             :    real(r8), intent(in) :: gcl(pcols,0:pver)    ! Cloud asymmetry parameter (liquid)
        2023             :    real(r8), intent(in) :: fcl(pcols,0:pver)    ! Cloud forward scattered fraction (liquid)
        2024             :    real(r8), intent(in) :: tauxci(pcols,0:pver) ! Cloud extinction optical depth (ice)
        2025             :    real(r8), intent(in) :: wci(pcols,0:pver)    ! Cloud single scattering albedo (ice)
        2026             :    real(r8), intent(in) :: gci(pcols,0:pver)    ! Cloud asymmetry parameter (ice)
        2027             :    real(r8), intent(in) :: fci(pcols,0:pver)    ! Cloud forward scattered fraction (ice)
        2028             :    real(r8), intent(inout) :: aer_tau(pcols,0:pver) ! Aerosol extinction optical depth
        2029             :    real(r8), intent(inout) :: aer_tau_w(pcols,0:pver)     ! Aerosol single scattering albedo * tau
        2030             :    real(r8), intent(inout) :: aer_tau_w_g(pcols,0:pver)     ! Aerosol asymmetry parameter * w * t
        2031             :    real(r8), intent(inout) :: aer_tau_w_f(pcols,0:pver)     ! Aerosol forward scattered fraction * w * tau
        2032             :
        2033             :    integer, intent(in) :: ndayc                 ! Number of daylight columns
        2034             :    integer, intent(in) :: ns                    ! Index of spectral interval
        2035             : !
        2036             : ! Input/Output arguments
        2037             : !
        2038             : ! Following variables are defined for each layer; 0 refers to extra
        2039             : ! layer above top of model:
        2040             : !
        2041             :    real(r8), intent(inout) :: rdir(nswbands,pcols,0:pver)   ! Layer reflectivity to direct rad
        2042             :    real(r8), intent(inout) :: rdif(nswbands,pcols,0:pver)   ! Layer reflectivity to diffuse rad
        2043             :    real(r8), intent(inout) :: tdir(nswbands,pcols,0:pver)   ! Layer transmission to direct rad
        2044             :    real(r8), intent(inout) :: tdif(nswbands,pcols,0:pver)   ! Layer transmission to diffuse rad
        2045             :    real(r8), intent(inout) :: explay(nswbands,pcols,0:pver) ! Solar beam exp transm for layer
        2046             : !
        2047             : ! Corresponding quantities for clear-skies
        2048             : !
        2049             :    real(r8), intent(inout) :: rdirc(nswbands,pcols,0:pver)  ! Clear layer reflec. to direct rad
        2050             :    real(r8), intent(inout) :: rdifc(nswbands,pcols,0:pver)  ! Clear layer reflec. to diffuse rad
        2051             :    real(r8), intent(inout) :: tdirc(nswbands,pcols,0:pver)  ! Clear layer trans. to direct rad
        2052             :    real(r8), intent(inout) :: tdifc(nswbands,pcols,0:pver)  ! Clear layer trans. to diffuse rad
        2053             :    real(r8), intent(inout) :: explayc(nswbands,pcols,0:pver)! Solar beam exp transm clear layer
        2054             : !
        2055             : !---------------------------Local variables-----------------------------
        2056             : !
        2057             :    integer i                 ! Column indices
        2058             :    integer k                 ! Level index
        2059             :    integer nn                ! Index of column loops (max=ndayc)
        2060             :
        2061             :    real(r8) taugab(pcols)        ! Layer total gas absorption optical depth
        2062             :    real(r8) tauray(pcols)        ! Layer rayleigh optical depth
        2063             :    real(r8) taucsc               ! Layer cloud scattering optical depth
        2064             :    real(r8) tautot               ! Total layer optical depth
        2065             :    real(r8) wtot                 ! Total layer single scatter albedo
        2066             :    real(r8) gtot                 ! Total layer asymmetry parameter
        2067             :    real(r8) ftot                 ! Total layer forward scatter fraction
        2068             :    real(r8) wtau                 !  rayleigh layer scattering optical depth
        2069             :    real(r8) wt                   !  layer total single scattering albedo
        2070             :    real(r8) ts                   !  layer scaled extinction optical depth
        2071             :    real(r8) ws                   !  layer scaled single scattering albedo
        2072             :    real(r8) gs                   !  layer scaled asymmetry parameter
        2073             : !
        2074             : !---------------------------Statement functions-------------------------
        2075             : !
        2076             : ! Statement functions and other local variables
        2077             : !
        2078             :    real(r8) alpha                ! Term in direct reflect and transmissivity
        2079             :    real(r8) gamma                ! Term in direct reflect and transmissivity
        2080             :    real(r8) el                   ! Term in alpha,gamma,n,u
        2081             :    real(r8) taus                 ! Scaled extinction optical depth
        2082             :    real(r8) omgs                 ! Scaled single particle scattering albedo
        2083             :    real(r8) asys                 ! Scaled asymmetry parameter
        2084             :    real(r8) u                    ! Term in diffuse reflect and
        2085             : !    transmissivity
        2086             :    real(r8) n                    ! Term in diffuse reflect and
        2087             : !    transmissivity
        2088             :    real(r8) lm                   ! Temporary for el
        2089             :    real(r8) ne                   ! Temporary for n
        2090             :    real(r8) w                    ! Dummy argument for statement function
        2091             :    real(r8) uu                   ! Dummy argument for statement function
        2092             :    real(r8) g                    ! Dummy argument for statement function
        2093             :    real(r8) e                    ! Dummy argument for statement function
        2094             :    real(r8) f                    ! Dummy argument for statement function
        2095             :    real(r8) t                    ! Dummy argument for statement function
        2096             :    real(r8) et                   ! Dummy argument for statement function
        2097             : !
        2098             : ! Intermediate terms for delta-eddington solution
        2099             : !
        2100             :    real(r8) alp                  ! Temporary for alpha
        2101             :    real(r8) gam                  ! Temporary for gamma
        2102             :    real(r8) ue                   ! Temporary for u
        2103             :    real(r8) arg                  ! Exponential argument
        2104             :    real(r8) extins               ! Extinction
        2105             :    real(r8) amg                  ! Alp - gam
        2106             :    real(r8) apg                  ! Alp + gam
        2107             : !
        2108             : ! ssa <=1 limit for aerosol
        2109             : !
        2110             :    real(r8) :: w_limited(pcols,0:pver)       ! Aerosol ssa (limited to < 0.999999)
        2111             :    real(r8) :: aer_g_limit(pcols,0:pver)     ! Aerosol tau_w_g (limited ssa)
        2112             :    real(r8) :: aer_f_limit(pcols,0:pver)     ! Aerosol tau_w_f (limited ssa)
        2113             : !
        2114             :    alpha(w,uu,g,e) = .75_r8*w*uu*((1._r8 + g*(1._r8-w))/(1._r8 - e*e*uu*uu))
        2115             :    gamma(w,uu,g,e) = .50_r8*w*((3._r8*g*(1._r8-w)*uu*uu + 1._r8)/(1._r8-e*e*uu*uu))
        2116             :    el(w,g)         = sqrt(3._r8*(1._r8-w)*(1._r8 - w*g))
        2117             :    taus(w,f,t)     = (1._r8 - w*f)*t
        2118             :    omgs(w,f)       = (1._r8 - f)*w/(1._r8 - w*f)
        2119             :    asys(g,f)       = (g - f)/(1._r8 - f)
        2120             :    u(w,g,e)        = 1.5_r8*(1._r8 - w*g)/e
        2121             :    n(uu,et)        = ((uu+1._r8)*(uu+1._r8)/et ) - ((uu-1._r8)*(uu-1._r8)*et)
        2122             : !
        2123             : !-----------------------------------------------------------------------
        2124             : !
        2125             : ! Compute layer radiative properties
        2126             : !
        2127             : ! Compute radiative properties (reflectivity and transmissivity for
        2128             : !    direct and diffuse radiation incident from above, under clear
        2129             : !    and cloudy conditions) and transmission of direct radiation
        2130             : !    (under clear and cloudy conditions) for each layer.
        2131             : !
        2132   208428556 :    do k=0,pver
        2133  3225176356 :       do i=1,ndayc
        2134  3016747800 :          if(aer_tau(i,k) > 0._r8) then !where(aer_tau > 0._r8)
        2135           0 :             aer_g_limit(i,k) = aer_tau_w_g(i,k) / aer_tau_w(i,k)
        2136           0 :             aer_f_limit(i,k) = aer_tau_w_f(i,k) / aer_tau_w(i,k)
        2137           0 :             aer_tau_w(i,k) = aer_tau(i,k) * min(aer_tau_w(i,k)/aer_tau(i,k) , 0.999999_r8)
        2138             :          else
        2139  3016747800 :             aer_tau(i,k) = 0._r8
        2140  3016747800 :             aer_tau_w(i,k) = 0._r8
        2141  3016747800 :             aer_g_limit(i,k) = 0._r8
        2142  3016747800 :             aer_f_limit(i,k) = 0._r8
        2143             :          endif
        2144  3016747800 :          aer_tau_w_g(i,k) = aer_tau_w(i,k) * aer_g_limit(i,k)
        2145  3217732479 :          aer_tau_w_f(i,k) = aer_tau_w(i,k) * aer_f_limit(i,k)
        2146             :       enddo
        2147             :    enddo
        2148             :
        2149   208428556 :    do k=0,pver
        2150  3225176356 :       do i=1,ndayc
        2151  3016747800 :             tauray(i) = trayoslp*(pflx(i,k+1)-pflx(i,k))
        2152  3016747800 :             taugab(i) = abh2o*uh2o(i,k) + abo3*uo3(i,k) + abco2*uco2(i,k) + abo2*uo2(i,k)
        2153  3016747800 :             tautot = tauxcl(i,k) + tauxci(i,k) + tauray(i) + taugab(i) + aer_tau(i,k)
        2154  3016747800 :             taucsc = tauxcl(i,k)*wcl(i,k) + tauxci(i,k)*wci(i,k) + aer_tau_w(i,k)
        2155  3016747800 :             wtau   = wray*tauray(i)
        2156  3016747800 :             wt     = wtau + taucsc
        2157  3016747800 :             wtot   = wt/tautot
        2158             :             gtot   = (wtau*gray + gcl(i,k)*wcl(i,k)*tauxcl(i,k) &
        2159  3016747800 :                      + gci(i,k)*wci(i,k)*tauxci(i,k) + aer_tau_w_g(i,k))/wt
        2160             :             ftot   = (wtau*fray + fcl(i,k)*wcl(i,k)*tauxcl(i,k) &
        2161  3016747800 :                      + fci(i,k)*wci(i,k)*tauxci(i,k) + aer_tau_w_f(i,k))/wt
        2162  3016747800 :             ts   = taus(wtot,ftot,tautot)
        2163  3016747800 :             ws   = omgs(wtot,ftot)
        2164  3016747800 :             gs   = asys(gtot,ftot)
        2165  3016747800 :             lm   = el(ws,gs)
        2166  3016747800 :             alp  = alpha(ws,coszrs(i),gs,lm)
        2167  3016747800 :             gam  = gamma(ws,coszrs(i),gs,lm)
        2168  3016747800 :             ue   = u(ws,gs,lm)
        2169             : !
        2170             : !     Limit argument of exponential to 25, in case lm very large:
        2171             : !
        2172  3016747800 :             arg  = min(lm*ts,25._r8)
        2173  3016747800 :             extins = exp(-arg)
        2174  3016747800 :             ne = n(ue,extins)
        2175  3016747800 :             rdif(ns,i,k) = (ue+1._r8)*(ue-1._r8)*(1._r8/extins - extins)/ne
        2176  3016747800 :             tdif(ns,i,k)   =   4._r8*ue/ne
        2177             : !
        2178             : !     Limit argument of exponential to 25, in case coszrs is very small:
        2179             : !
        2180  3016747800 :             arg       = min(ts/coszrs(i),25._r8)
        2181  3016747800 :             explay(ns,i,k) = exp(-arg)
        2182  3016747800 :             apg = alp + gam
        2183  3016747800 :             amg = alp - gam
        2184  3016747800 :             rdir(ns,i,k) = amg*(tdif(ns,i,k)*explay(ns,i,k)-1._r8) + apg*rdif(ns,i,k)
        2185  3016747800 :             tdir(ns,i,k) = apg*tdif(ns,i,k) + (amg*rdif(ns,i,k)-(apg-1._r8))*explay(ns,i,k)
        2186             : !
        2187             : !     Under rare conditions, reflectivies and transmissivities can be
        2188             : !     negative; zero out any negative values
        2189             : !
        2190  3016747800 :             rdir(ns,i,k) = max(rdir(ns,i,k),0.0_r8)
        2191  3016747800 :             tdir(ns,i,k) = max(tdir(ns,i,k),0.0_r8)
        2192  3016747800 :             rdif(ns,i,k) = max(rdif(ns,i,k),0.0_r8)
        2193  3016747800 :             tdif(ns,i,k) = max(tdif(ns,i,k),0.0_r8)
        2194             : !
        2195             : !     Clear-sky calculation
        2196             : !
        2197  3217732479 :             if (tauxcl(i,k) == 0.0_r8 .and. tauxci(i,k) == 0.0_r8) then
        2198             :
        2199  2281020946 :                rdirc(ns,i,k) = rdir(ns,i,k)
        2200  2281020946 :                tdirc(ns,i,k) = tdir(ns,i,k)
        2201  2281020946 :                rdifc(ns,i,k) = rdif(ns,i,k)
        2202  2281020946 :                tdifc(ns,i,k) = tdif(ns,i,k)
        2203  2281020946 :                explayc(ns,i,k) = explay(ns,i,k)
        2204             :             else
        2205   735726854 :                tautot = tauray(i) + taugab(i) + aer_tau(i,k)
        2206   735726854 :                taucsc = aer_tau_w(i,k)
        2207             : !
        2208             : ! wtau already computed for all-sky
        2209             : !
        2210   735726854 :                wt     = wtau + taucsc
        2211   735726854 :                wtot   = wt/tautot
        2212   735726854 :                gtot   = (wtau*gray + aer_tau_w_g(i,k))/wt
        2213   735726854 :                ftot   = (wtau*fray + aer_tau_w_f(i,k))/wt
        2214   735726854 :                ts   = taus(wtot,ftot,tautot)
        2215   735726854 :                ws   = omgs(wtot,ftot)
        2216   735726854 :                gs   = asys(gtot,ftot)
        2217   735726854 :                lm   = el(ws,gs)
        2218   735726854 :                alp  = alpha(ws,coszrs(i),gs,lm)
        2219   735726854 :                gam  = gamma(ws,coszrs(i),gs,lm)
        2220   735726854 :                ue   = u(ws,gs,lm)
        2221             : !
        2222             : !     Limit argument of exponential to 25, in case lm very large:
        2223             : !
        2224   735726854 :                arg  = min(lm*ts,25._r8)
        2225   735726854 :                extins = exp(-arg)
        2226   735726854 :                ne = n(ue,extins)
        2227   735726854 :                rdifc(ns,i,k) = (ue+1._r8)*(ue-1._r8)*(1._r8/extins - extins)/ne
        2228   735726854 :                tdifc(ns,i,k)   =   4._r8*ue/ne
        2229             : !
        2230             : !     Limit argument of exponential to 25, in case coszrs is very small:
        2231             : !
        2232   735726854 :                arg       = min(ts/coszrs(i),25._r8)
        2233   735726854 :                explayc(ns,i,k) = exp(-arg)
        2234   735726854 :                apg = alp + gam
        2235   735726854 :                amg = alp - gam
        2236             :                rdirc(ns,i,k) = amg*(tdifc(ns,i,k)*explayc(ns,i,k)-1._r8)+ &
        2237   735726854 :                                apg*rdifc(ns,i,k)
        2238             :                tdirc(ns,i,k) = apg*tdifc(ns,i,k) + (amg*rdifc(ns,i,k) - (apg-1._r8))* &
        2239   735726854 :                                explayc(ns,i,k)
        2240             : !
        2241             : !     Under rare conditions, reflectivies and transmissivities can be
        2242             : !     negative; zero out any negative values
        2243             : !
        2244   735726854 :                rdirc(ns,i,k) = max(rdirc(ns,i,k),0.0_r8)
        2245   735726854 :                tdirc(ns,i,k) = max(tdirc(ns,i,k),0.0_r8)
        2246   735726854 :                rdifc(ns,i,k) = max(rdifc(ns,i,k),0.0_r8)
        2247   735726854 :                tdifc(ns,i,k) = max(tdifc(ns,i,k),0.0_r8)
        2248             :             end if
        2249             :          end do
        2250             :    end do
        2251             :
        2252      749232 : end subroutine raddedmx
        2253             :
        2254             : !-------------------------------------------------------------------------------
        2255             :
        2256        1536 : subroutine radsw_init(gravx)
        2257             : !-----------------------------------------------------------------------
        2258             : !
        2259             : ! Purpose:
        2260             : ! Initialize various constants for radiation scheme; note that
        2261             : ! the radiation scheme uses cgs units.
        2262             : !
        2263             : ! Author: W. Collins (H2O parameterization) and J. Kiehl
        2264             : !
        2265             : !-----------------------------------------------------------------------
        2266             : !
        2267             : ! Input arguments
        2268             : !
        2269             :    real(r8), intent(in) :: gravx      ! Acceleration of gravity (MKS)
        2270             :
        2271             :    real(r8), parameter :: ref_tsi = 1367._r8 ! value supplied by Dan Marsh -- see solvar_woods.F90
        2272             : !
        2273             : !-----------------------------------------------------------------------
        2274             : !
        2275             : ! Set general radiation consts; convert to cgs units where appropriate:
        2276             : !
        2277        1536 :    gravit  =  100._r8*gravx
        2278        1536 :    rga     =  1._r8/gravit
        2279        1536 :    sslp    =  1.013250e6_r8
        2280             :
        2281        1536 : end subroutine radsw_init
        2282             :
        2283             : !-------------------------------------------------------------------------------
        2284             :
        2285             : end module radsw

Generated by: LCOV version 1.14