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

              Line data    Source code
           1             :
           2             : module radlw
           3             : !-----------------------------------------------------------------------
           4             : !
           5             : ! Purpose: Longwave 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 scamMod,           only: single_column, scm_crm_mode
          12             : use radconstants,      only: nlwbands
          13             : implicit none
          14             :
          15             : private
          16             : save
          17             :
          18             : ! Public methods
          19             :
          20             : public ::&
          21             :    radlw_init,   &! initialize constants
          22             :    radclwmx       ! driver for longwave radiation code
          23             :
          24             : ! Private module data
          25             :
          26             : real(r8) :: gravit_cgs     ! Acceleration of gravity (cm/s**2)
          27             : real(r8) :: stebol_cgs     ! Stefan-Boltzmann constant (CGS)
          28             :
          29             : !===============================================================================
          30             : CONTAINS
          31             : !===============================================================================
          32             :
          33      749232 : subroutine radclwmx(lchnk   ,ncol    ,doabsems                  , &
          34             :                     lwupcgs ,tnm     ,qnm     ,o3      , &
          35             :                     pmid    ,pint    ,pmln    ,piln    ,          &
          36             :                              n2o     ,ch4     ,cfc11   ,cfc12   , &
          37             :                     cld     ,emis    ,pmxrgn  ,nmxrgn  ,qrl,qrlc, &
          38             :                     flns    ,flnt    ,flnsc   ,flntc   ,flwds   , &
          39             :                     fldsc   ,flut    ,flutc   , &
          40             :                     fnl     ,fcnl    ,co2mmr,  odap_aer)
          41             : !-----------------------------------------------------------------------
          42             : !
          43             : ! Purpose:
          44             : ! Compute longwave radiation heating rates and boundary fluxes
          45             : !
          46             : ! Method:
          47             : ! Uses broad band absorptivity/emissivity method to compute clear sky;
          48             : ! assumes randomly overlapped clouds with variable cloud emissivity to
          49             : ! include effects of clouds.
          50             : !
          51             : ! Computes clear sky absorptivity/emissivity at lower frequency (in
          52             : ! general) than the model radiation frequency; uses previously computed
          53             : ! and stored values for efficiency
          54             : !
          55             : ! Note: This subroutine contains vertical indexing which proceeds
          56             : !       from bottom to top rather than the top to bottom indexing
          57             : !       used in the rest of the model.
          58             : !
          59             : ! Author: B. Collins
          60             : !
          61             : !-----------------------------------------------------------------------
          62             :    use radae,           only: nbands, radems, radabs, radtpl, abstot_3d, &
          63             :                               absnxt_3d, emstot_3d, ntoplw, radoz2, trcpth
          64             :    use cam_history,     only: outfld
          65             :    use quicksort,       only: quick_sort
          66             :    use radconstants,    only: nlwbands
          67             :    use phys_control,    only: cam_physpkg_is
          68             :    use ref_pres,        only: trop_cloud_top_lev
          69             :
          70             :    integer,  parameter :: pverp2=pver+2, pverp3=pver+3, pverp4=pver+4
          71             :    real(r8), parameter :: cldmin = 1.0e-80_r8
          72             :
          73             : !------------------------------Arguments--------------------------------
          74             : !
          75             : ! Input arguments
          76             : !
          77             :    integer, intent(in) :: lchnk                 ! chunk identifier
          78             :    integer, intent(in) :: ncol                  ! number of atmospheric columns
          79             :    logical, intent(in) :: doabsems              ! True => abs/emiss calculation this timestep
          80             :
          81             : !    maximally overlapped region.
          82             : !    0->pmxrgn(i,1) is range of pmid for
          83             : !    1st region, pmxrgn(i,1)->pmxrgn(i,2) for
          84             : !    2nd region, etc
          85             :    integer, intent(in) :: nmxrgn(pcols)         ! Number of maximally overlapped regions
          86             :
          87             :    real(r8), intent(in) :: pmxrgn(pcols,pverp)  ! Maximum values of pmid for each
          88             :    real(r8), intent(in) :: lwupcgs(pcols)       ! Longwave up flux in CGS units
          89             : !
          90             : ! Input arguments which are only passed to other routines
          91             : !
          92             :    real(r8), intent(in) :: tnm(pcols,pver)      ! Level temperature
          93             :    real(r8), intent(in) :: qnm(pcols,pver)      ! Level moisture field
          94             :    real(r8), intent(in) :: o3(pcols,pver)       ! ozone mass mixing ratio
          95             :    real(r8), intent(in) :: pmid(pcols,pver)     ! Level pressure
          96             :    real(r8), intent(in) :: pint(pcols,pverp)    ! Model interface pressure
          97             :    real(r8), intent(in) :: pmln(pcols,pver)     ! Ln(pmid)
          98             :    real(r8), intent(in) :: piln(pcols,pverp)    ! Ln(pint)
          99             :    real(r8), intent(in) :: co2mmr(pcols)        ! co2 column mean mass mixing ratio
         100             :    real(r8), intent(in) :: n2o(pcols,pver)      ! nitrous oxide mass mixing ratio
         101             :    real(r8), intent(in) :: ch4(pcols,pver)      ! methane mass mixing ratio
         102             :    real(r8), intent(in) :: cfc11(pcols,pver)    ! cfc11 mass mixing ratio
         103             :    real(r8), intent(in) :: cfc12(pcols,pver)    ! cfc12 mass mixing ratio
         104             :    real(r8), intent(in) :: cld(pcols,pver)      ! Cloud cover
         105             :    real(r8), intent(in) :: emis(pcols,pver)     ! Cloud emissivity
         106             :
         107             : ! [fraction] absorbtion optical depth, cumulative from top
         108             :   real(r8), intent(in) ::  odap_aer(pcols,pver,nlwbands)
         109             :
         110             :
         111             : !
         112             : ! Output arguments
         113             : !
         114             :    real(r8), intent(out) :: qrl (pcols,pver)     ! Longwave heating rate
         115             :    real(r8), intent(out) :: qrlc(pcols,pver)     ! Clearsky longwave heating rate
         116             :    real(r8), intent(out) :: flns(pcols)          ! Surface cooling flux
         117             :    real(r8), intent(out) :: flnt(pcols)          ! Net outgoing flux
         118             :    real(r8), intent(out) :: flut(pcols)          ! Upward flux at top of model
         119             :    real(r8), intent(out) :: flnsc(pcols)         ! Clear sky surface cooing
         120             :    real(r8), intent(out) :: flntc(pcols)         ! Net clear sky outgoing flux
         121             :    real(r8), intent(out) :: flutc(pcols)         ! Upward clear-sky flux at top of model
         122             :    real(r8), intent(out) :: flwds(pcols)         ! Down longwave flux at surface
         123             :    real(r8), intent(out) :: fldsc(pcols)         ! Down clear-sky longwave flux at surface
         124             :    real(r8),intent(out)  :: fcnl(pcols,pverp)    ! clear sky net flux at interfaces
         125             :    real(r8),intent(out)  :: fnl(pcols,pverp)     ! net flux at interfaces
         126             : !
         127             : !---------------------------Local variables-----------------------------
         128             : !
         129             :    ! Implicit save here is fine since this should have the same value for
         130             :    ! the entire run.
         131             :    integer :: ntopcld = 2    ! Lowest layer without clouds.
         132             :                              ! Shouldn't this be turned off by default?
         133             :
         134             :    integer i                 ! Longitude index
         135             :    integer ilon              ! Longitude index
         136             :    integer ii                ! Longitude index
         137             :    integer iimx              ! Longitude index (max overlap)
         138             :    integer k                 ! Level index
         139             :    integer k1                ! Level index
         140             :    integer k2                ! Level index
         141             :    integer k3                ! Level index
         142             :    integer km                ! Level index
         143             :    integer km1               ! Level index
         144             :    integer km3               ! Level index
         145             :    integer km4               ! Level index
         146             :    integer irgn              ! Index for max-overlap regions
         147             :    integer l                 ! Index for clouds to overlap
         148             :    integer l1                ! Index for clouds to overlap
         149             :    integer n                 ! Counter
         150             :
         151             : !
         152             :    real(r8) :: plco2(pcols,pverp)   ! Path length co2
         153             :    real(r8) :: plh2o(pcols,pverp)   ! Path length h2o
         154             :    real(r8) tmp(pcols)           ! Temporary workspace
         155             :    real(r8) tmp2(pcols)          ! Temporary workspace
         156             :    real(r8) tmp3(0:pverp)        ! Temporary workspace
         157             :    real(r8) tmp4                 ! Temporary workspace
         158             :    real(r8) tfdl                 ! Temporary workspace
         159             :    real(r8) tful                 ! Temporary workspace
         160             :    real(r8) absbt(pcols)         ! Downward emission at model top
         161             :    real(r8) plol(pcols,pverp)    ! O3 pressure wghted path length
         162             :    real(r8) plos(pcols,pverp)    ! O3 path length
         163             :    real(r8) co2em(pcols,pverp)   ! Layer co2 normalized planck funct. derivative
         164             :    real(r8) co2eml(pcols,pver)   ! Interface co2 normalized planck funct. deriv.
         165             :    real(r8) delt(pcols)          ! Diff t**4 mid layer to top interface
         166             :    real(r8) delt1(pcols)         ! Diff t**4 lower intrfc to mid layer
         167             :    real(r8) bk1(pcols)           ! Absrptvty for vertical quadrature
         168             :    real(r8) bk2(pcols)           ! Absrptvty for vertical quadrature
         169             :    real(r8) cldp(pcols,pverp)    ! Cloud cover with extra layer
         170             :    real(r8) ful(pcols,pverp)     ! Total upwards longwave flux
         171             :    real(r8) fsul(pcols,pverp)    ! Clear sky upwards longwave flux
         172             :    real(r8) fdl(pcols,pverp)     ! Total downwards longwave flux
         173             :    real(r8) fsdl(pcols,pverp)    ! Clear sky downwards longwv flux
         174             :    real(r8) fclb4(pcols,-1:pver)    ! Sig t**4 for cld bottom interfc
         175             :    real(r8) fclt4(pcols,0:pver)    ! Sig t**4 for cloud top interfc
         176             :    real(r8) s(pcols,pverp,pverp) ! Flx integral sum
         177             :    real(r8) tplnka(pcols,pverp)  ! Planck fnctn temperature
         178             :    real(r8) s2c(pcols,pverp)     ! H2o cont amount
         179             :    real(r8) tcg(pcols,pverp)     ! H2o-mass-wgted temp. (Curtis-Godson approx.)
         180             :    real(r8) w(pcols,pverp)       ! H2o path
         181             :    real(r8) tplnke(pcols)        ! Planck fnctn temperature
         182             :    real(r8) h2otr(pcols,pverp)   ! H2o trnmsn for o3 overlap
         183             :    real(r8) co2t(pcols,pverp)    ! Prs wghted temperature path
         184             :    real(r8) tint(pcols,pverp)    ! Interface temperature
         185             :    real(r8) tint4(pcols,pverp)   ! Interface temperature**4
         186             :    real(r8) tlayr(pcols,pverp)   ! Level temperature
         187             :    real(r8) tlayr4(pcols,pverp)  ! Level temperature**4
         188             :    real(r8) plh2ob(nbands,pcols,pverp)! Pressure weighted h2o path with
         189             :                                       !    Hulst-Curtis-Godson temp. factor
         190             :                                       !    for H2O bands
         191             :    real(r8) wb(nbands,pcols,pverp)    ! H2o path length with
         192             :                                       !    Hulst-Curtis-Godson temp. factor
         193             :                                       !    for H2O bands
         194             :
         195             :    real(r8) cld0                 ! previous cloud amt (for max overlap)
         196             :    real(r8) cld1                 ! next cloud amt (for max overlap)
         197             :    real(r8) emx(0:pverp)         ! Emissivity factors (max overlap)
         198             :    real(r8) emx0                 ! Emissivity factors for BCs (max overlap)
         199             :    real(r8) trans                ! 1 - emis
         200             :    real(r8) asort(pver)          ! 1 - cloud amounts to be sorted for max ovrlp.
         201             :    real(r8) atmp                 ! Temporary storage for sort when nxs = 2
         202             :    real(r8) maxcld(pcols)        ! Maximum cloud at any layer
         203             :
         204             :    integer indx(pcols)       ! index vector of gathered array values
         205             :    integer indxmx(pcols+1,pverp)! index vector of gathered array values
         206             : !   integer indxmx(pcols,pverp)! index vector of gathered array values
         207             : !    (max overlap)
         208             :    integer nrgn(pcols)       ! Number of max overlap regions at longitude
         209             :    integer npts              ! number of values satisfying some criterion
         210             :    integer ncolmx(pverp)     ! number of columns with clds in region
         211             :    integer kx1(pcols,pverp)  ! Level index for top of max-overlap region
         212             :    integer kx2(pcols,0:pverp)! Level index for bottom of max-overlap region
         213             :    integer kxs(0:pverp,pcols,pverp)! Level indices for cld layers sorted by cld()
         214             : !    in descending order
         215             :    integer nxs(pcols,pverp)  ! Number of cloudy layers between kx1 and kx2
         216             :    integer nxsk              ! Number of cloudy layers between (kx1/kx2)&k
         217             :    integer ksort(0:pverp)    ! Level indices of cloud amounts to be sorted
         218             : !    for max ovrlp. calculation
         219             :    integer ktmp              ! Temporary storage for sort when nxs = 2
         220             :
         221             : !
         222             : ! Pointer variables to 3d structures
         223             : !
         224      749232 :    real(r8), pointer :: abstot(:,:,:)
         225      749232 :    real(r8), pointer :: absnxt(:,:,:)
         226      749232 :    real(r8), pointer :: emstot(:,:)
         227             :
         228             : ! [fraction] Total transmission between interfaces k1 and k2
         229             :    real(r8) ::  aer_trn_ttl(pcols,pverp,pverp,nlwbands)
         230             : !
         231             : ! Trace gas variables
         232             : !
         233             :    real(r8) ucfc11(pcols,pverp)  ! CFC11 path length
         234             :    real(r8) ucfc12(pcols,pverp)  ! CFC12 path length
         235             :    real(r8) un2o0(pcols,pverp)   ! N2O path length
         236             :    real(r8) un2o1(pcols,pverp)   ! N2O path length (hot band)
         237             :    real(r8) uch4(pcols,pverp)    ! CH4 path length
         238             :    real(r8) uco211(pcols,pverp)  ! CO2 9.4 micron band path length
         239             :    real(r8) uco212(pcols,pverp)  ! CO2 9.4 micron band path length
         240             :    real(r8) uco213(pcols,pverp)  ! CO2 9.4 micron band path length
         241             :    real(r8) uco221(pcols,pverp)  ! CO2 10.4 micron band path length
         242             :    real(r8) uco222(pcols,pverp)  ! CO2 10.4 micron band path length
         243             :    real(r8) uco223(pcols,pverp)  ! CO2 10.4 micron band path length
         244             :    real(r8) bn2o0(pcols,pverp)   ! pressure factor for n2o
         245             :    real(r8) bn2o1(pcols,pverp)   ! pressure factor for n2o
         246             :    real(r8) bch4(pcols,pverp)    ! pressure factor for ch4
         247             :    real(r8) uptype(pcols,pverp)  ! p-type continuum path length
         248             :    real(r8) abplnk1(14,pcols,pverp)  ! non-nearest layer Plack factor
         249             :    real(r8) abplnk2(14,pcols,pverp)  ! nearest layer factor
         250             : !
         251             : !
         252             : !-----------------------------------------------------------------------
         253             : !
         254             : !
         255             : ! Set pointer variables
         256             : !
         257      749232 :    abstot => abstot_3d(:,:,:,lchnk)
         258      749232 :    absnxt => absnxt_3d(:,:,:,lchnk)
         259      749232 :    emstot => emstot_3d(:,:,lchnk)
         260             :
         261             : !
         262             : ! Calculate some temperatures needed to derive absorptivity and
         263             : ! emissivity, as well as some h2o path lengths
         264             : !
         265             :    call radtpl(ncol    ,                                     &
         266             :                tnm     ,lwupcgs ,qnm     ,pint    ,plco2   ,plh2o   , &
         267             :                tplnka  ,s2c     ,tcg     ,w       ,tplnke  , &
         268             :                tint    ,tint4   ,tlayr   ,tlayr4  ,pmln    , &
         269      749232 :                piln    ,plh2ob  ,wb      ,co2mmr)
         270             :
         271             :
         272      749232 :    if (doabsems) then
         273             : !
         274             : ! Compute ozone path lengths at frequency of a/e calculation.
         275             : !
         276       68112 :       call radoz2(ncol, o3, pint, plol, plos)
         277             : !
         278             : ! Compute trace gas path lengths
         279             : !
         280             :       call trcpth(ncol                                        , &
         281             :                   tnm     ,pint    ,cfc11   ,cfc12   ,n2o     , &
         282             :                   ch4     ,qnm     ,ucfc11  ,ucfc12  ,un2o0   , &
         283             :                   un2o1   ,uch4    ,uco211  ,uco212  ,uco213  , &
         284             :                   uco221  ,uco222  ,uco223  ,bn2o0   ,bn2o1   , &
         285       68112 :                   bch4    ,uptype  ,co2mmr)
         286             : !
         287             : ! Compute transmission through aerosols from (absorption) optical depths
         288             : !
         289       68112 :       call aer_trans_from_od(ncol, odap_aer, aer_trn_ttl)
         290             : !
         291             : ! Compute total emissivity:
         292             : !
         293             :       call radems(lchnk   ,ncol    ,                            &
         294             :                   s2c     ,tcg     ,w       ,tplnke  ,plh2o   , &
         295             :                   pint    ,plco2   ,tint    ,tint4   ,tlayr   , &
         296             :                   tlayr4  ,plol    ,plos    ,ucfc11  ,ucfc12  , &
         297             :                   un2o0   ,un2o1   ,uch4    ,uco211  ,uco212  , &
         298             :                   uco213  ,uco221  ,uco222  ,uco223  ,uptype  , &
         299             :                   bn2o0   ,bn2o1   ,bch4    ,co2em   ,co2eml  , &
         300             :                   co2t    ,h2otr   ,abplnk1 ,abplnk2 ,emstot  , &
         301             :                   plh2ob  ,wb      , &
         302       68112 :                   aer_trn_ttl, co2mmr)
         303             : !
         304             : ! Compute total absorptivity:
         305             : !
         306             :
         307             :       call radabs(lchnk   ,ncol    ,                            &
         308             :                   pmid    ,pint    ,co2em   ,co2eml  ,tplnka  , &
         309             :                   s2c     ,tcg     ,w       ,h2otr   ,plco2   , &
         310             :                   plh2o   ,co2t    ,tint    ,tlayr   ,plol    , &
         311             :                   plos    ,pmln    ,piln    ,ucfc11  ,ucfc12  , &
         312             :                   un2o0   ,un2o1   ,uch4    ,uco211  ,uco212  , &
         313             :                   uco213  ,uco221  ,uco222  ,uco223  ,uptype  , &
         314             :                   bn2o0   ,bn2o1   ,bch4    ,abplnk1 ,abplnk2 , &
         315             :                   abstot  ,absnxt  ,plh2ob  ,wb      , &
         316       68112 :                   odap_aer,aer_trn_ttl, co2mmr)
         317             :    end if
         318             : !
         319             : ! Compute sums used in integrals (all longitude points)
         320             : !
         321             : ! Definition of bk1 & bk2 depends on finite differencing.  for
         322             : ! trapezoidal rule bk1=bk2. trapezoidal rule applied for nonadjacent
         323             : ! layers only.
         324             : !
         325             : ! delt=t**4 in layer above current sigma level km.
         326             : ! delt1=t**4 in layer below current sigma level km.
         327             : !
         328    12510432 :    do i=1,ncol
         329    11761200 :       delt(i) = tint4(i,pver) - tlayr4(i,pverp)
         330    11761200 :       delt1(i) = tlayr4(i,pverp) - tint4(i,pverp)
         331    11761200 :       s(i,pverp,pverp) = stebol_cgs*(delt1(i)*absnxt(i,pver,1) + delt (i)*absnxt(i,pver,4))
         332    12510432 :       s(i,pver,pverp)  = stebol_cgs*(delt (i)*absnxt(i,pver,2) + delt1(i)*absnxt(i,pver,3))
         333             :    end do
         334    19480032 :    do k=ntoplw,pver-1
         335   313510032 :       do i=1,ncol
         336   294030000 :          bk2(i) = (abstot_3d(i,k,pver,lchnk) + abstot_3d(i,k,pverp,lchnk))*0.5_r8
         337   294030000 :          bk1(i) = bk2(i)
         338   312760800 :          s(i,k,pverp) = stebol_cgs*(bk2(i)*delt(i) + bk1(i)*delt1(i))
         339             :       end do
         340             :    end do
         341             : !
         342             : ! All k, km>1
         343             : !
         344    19480032 :    do km=pver,ntoplw+1,-1
         345   312760800 :       do i=1,ncol
         346   294030000 :          delt(i)  = tint4(i,km-1) - tlayr4(i,km)
         347   312760800 :          delt1(i) = tlayr4(i,km) - tint4(i,km)
         348             :       end do
         349             : !CSD$ PARALLEL DO PRIVATE( i, k, bk1, bk2 )
         350   525211632 :       do k=pverp,ntoplw,-1
         351   505731600 :          if (k == km) then
         352   312760800 :             do i=1,ncol
         353   294030000 :                bk2(i) = absnxt(i,km-1,4)
         354   312760800 :                bk1(i) = absnxt(i,km-1,1)
         355             :             end do
         356   487000800 :          else if (k == km-1) then
         357   312760800 :             do i=1,ncol
         358   294030000 :                bk2(i) = absnxt(i,km-1,2)
         359   312760800 :                bk1(i) = absnxt(i,km-1,3)
         360             :             end do
         361             :          else
         362  7819020000 :             do i=1,ncol
         363  7350750000 :                bk2(i) = (abstot_3d(i,k,km-1,lchnk) + abstot_3d(i,k,km,lchnk))*0.5_r8
         364  7819020000 :                bk1(i) = bk2(i)
         365             :             end do
         366             :          end if
         367  8463272400 :          do i=1,ncol
         368  8444541600 :             s(i,k,km) = s(i,k,km+1) + stebol_cgs*(bk2(i)*delt(i) + bk1(i)*delt1(i))
         369             :          end do
         370             :       end do
         371             : !CSD$ END PARALLEL DO
         372             :    end do
         373             : !
         374             : ! Computation of clear sky fluxes always set first level of fsul
         375             : !
         376    12510432 :    do i=1,ncol
         377    12510432 :       fsul(i,pverp) = lwupcgs(i)
         378             :    end do
         379             : !
         380             : ! Downward clear sky fluxes store intermediate quantities in down flux
         381             : ! Initialize fluxes to clear sky values.
         382             : !
         383    12510432 :    do i=1,ncol
         384    11761200 :       tmp(i) = fsul(i,pverp) - stebol_cgs*tint4(i,pverp)
         385    11761200 :       fsul(i,ntoplw) = fsul(i,pverp) - abstot_3d(i,ntoplw,pverp,lchnk)*tmp(i) + s(i,ntoplw,ntoplw+1)
         386    12510432 :       fsdl(i,ntoplw) = stebol_cgs*(tplnke(i)**4)*emstot(i,ntoplw)
         387             :    end do
         388             : !
         389             : ! fsdl(i,pverp) assumes isothermal layer
         390             : !
         391    19480032 :    do k=ntoplw+1,pver
         392   313510032 :       do i=1,ncol
         393   294030000 :          fsul(i,k) = fsul(i,pverp) - abstot_3d(i,k,pverp,lchnk)*tmp(i) + s(i,k,k+1)
         394   312760800 :          fsdl(i,k) = stebol_cgs*(tplnke(i)**4)*emstot(i,k) - (s(i,k,ntoplw+1) - s(i,k,k+1))
         395             :       end do
         396             :    end do
         397             : !
         398             : ! Store the downward emission from level 1 = total gas emission * sigma
         399             : ! t**4.  fsdl does not yet include all terms
         400             : !
         401    12510432 :    do i=1,ncol
         402    11761200 :       absbt(i) = stebol_cgs*(tplnke(i)**4)*emstot(i,pverp)
         403    12510432 :       fsdl(i,pverp) = absbt(i) - s(i,pverp,ntoplw+1)
         404             :    end do
         405             :
         406    20978496 :    do k = ntoplw,pverp
         407   338530896 :       do i = 1,ncol
         408   337781664 :          fcnl(i,k) = fsul(i,k) - fsdl(i,k)
         409             :       end do
         410             :    end do
         411             : !
         412             : !----------------------------------------------------------------------
         413             : ! Modifications for clouds -- max/random overlap assumption
         414             : !
         415             : ! The column is divided into sets of adjacent layers, called regions,
         416             : !   in which the clouds are maximally overlapped.  The clouds are
         417             : !   randomly overlapped between different regions.  The number of
         418             : !   regions in a column is set by nmxrgn, and the range of pressures
         419             : !   included in each region is set by pmxrgn.  The max/random overlap
         420             : !   can be written in terms of the solutions of random overlap with
         421             : !   cloud amounts = 1.  The random overlap assumption is equivalent to
         422             : !   setting the flux boundary conditions (BCs) at the edges of each region
         423             : !   equal to the mean all-sky flux at those boundaries.  Since the
         424             : !   emissivity array for propogating BCs is only computed for the
         425             : !   TOA BC, the flux BCs elsewhere in the atmosphere have to be formulated
         426             : !   in terms of solutions to the random overlap equations.  This is done
         427             : !   by writing the flux BCs as the sum of a clear-sky flux and emission
         428             : !   from a cloud outside the region weighted by an emissivity.  This
         429             : !   emissivity is determined from the location of the cloud and the
         430             : !   flux BC.
         431             : !
         432             : ! Copy cloud amounts to buffer with extra layer (needed for overlap logic)
         433             : !
         434             :
         435      749232 :    ntopcld = max(ntopcld, trop_cloud_top_lev)
         436             :
         437    25770096 :    cldp(:ncol,1:ntopcld) = 0.0_r8
         438   300999600 :    cldp(:ncol,ntopcld+1:pver) = cld(:ncol,ntopcld+1:pver)
         439    12510432 :    cldp(:ncol,pverp) = 0.0_r8
         440             : !
         441             : !
         442             : ! Select only those locations where there are no clouds
         443             : !    (maximum cloud fraction <= 1.e-3 treated as clear)
         444             : !    Set all-sky fluxes to clear-sky values.
         445             : !
         446      749232 :    maxcld(1:ncol) = maxval(cldp(1:ncol,ntoplw:pver),dim=2)
         447             :
         448      749232 :    npts = 0
         449    12510432 :    do i=1,ncol
         450    12510432 :       if (maxcld(i) < cldmin) then
         451      316626 :          npts = npts + 1
         452      316626 :          indx(npts) = i
         453             :       end if
         454             :    end do
         455             :
         456     1065858 :    do ii = 1, npts
         457      316626 :       i = indx(ii)
         458     9614760 :       do k = ntoplw, pverp
         459     8548902 :          fdl(i,k) = fsdl(i,k)
         460     8865528 :          ful(i,k) = fsul(i,k)
         461             :       end do
         462             :    end do
         463             : !
         464             : ! Select only those locations where there are clouds
         465             : !
         466             :    npts = 0
         467    12510432 :    do i=1,ncol
         468    12510432 :       if (maxcld(i) >= cldmin) then
         469    11444574 :          npts = npts + 1
         470    11444574 :          indx(npts) = i
         471             :       end if
         472             :    end do
         473             :
         474             : !
         475             : ! Initialize all-sky fluxes. fdl(i,1) & ful(i,pverp) are boundary conditions
         476             : !
         477    12193806 :    do ii = 1, npts
         478    11444574 :       i = indx(ii)
         479    11444574 :       fdl(i,ntoplw) = fsdl(i,ntoplw)
         480    11444574 :       fdl(i,pverp)  = 0.0_r8
         481    11444574 :       ful(i,ntoplw) = 0.0_r8
         482    11444574 :       ful(i,pverp)  = fsul(i,pverp)
         483   297558924 :       do k = ntoplw+1, pver
         484   286114350 :          fdl(i,k) = 0.0_r8
         485   297558924 :          ful(i,k) = 0.0_r8
         486             :       end do
         487             : !
         488             : ! Initialize Planck emission from layer boundaries
         489             : !
         490   309003498 :       do k = ntoplw, pver
         491   297558924 :          fclt4(i,k-1) = stebol_cgs*tint4(i,k)
         492   309003498 :          fclb4(i,k-1) = stebol_cgs*tint4(i,k+1)
         493             :       enddo
         494    11444574 :       fclb4(i,ntoplw-2) =  stebol_cgs*tint4(i,ntoplw)
         495    11444574 :       fclt4(i,pver)     = stebol_cgs*tint4(i,pverp)
         496             : !
         497             : ! Initialize indices for layers to be max-overlapped
         498             : !
         499    40753317 :       do irgn = 0, nmxrgn(i)
         500    40753317 :          kx2(i,irgn) = ntoplw-1
         501             :       end do
         502    12193806 :       nrgn(i) = 0
         503             :    end do
         504             :
         505             : !----------------------------------------------------------------------
         506             : ! INDEX CALCULATIONS FOR MAX OVERLAP
         507             :
         508             : !CSD$ PARALLEL DO PRIVATE( ii, ilon, irgn, n, k1, k2, k, indxmx, ncolmx, iimx, i, ksort ) &
         509             : !CSD$ PRIVATE( asort, ktmp, atmp, km1, km4, k3, emx0, nxsk, emx, cld0 ) &
         510             : !CSD$ PRIVATE( tmp4, l1, tmp3, tfdl, l, cld1, trans, km3, tful )
         511    12193806 :    do ii = 1, npts
         512    11444574 :       ilon = indx(ii)
         513             :
         514             : !
         515             : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
         516             : !
         517    29308743 :       do irgn = 1, nmxrgn(ilon)
         518             : !
         519             : ! Calculate min/max layer indices inside region.
         520             : !
         521    17864169 :          n = 0
         522    17864169 :          if (kx2(ilon,irgn-1) < pver) then
         523    17864169 :             nrgn(ilon) = irgn
         524    17864169 :             k1 = kx2(ilon,irgn-1)+1
         525    17864169 :             kx1(ilon,irgn) = k1
         526    17864169 :             kx2(ilon,irgn) = k1 - 1
         527    89740248 :             do k2 = pver, k1, -1
         528    89740248 :                if (pmid(ilon,k2) <= pmxrgn(ilon,irgn)) then
         529    17864169 :                   kx2(ilon,irgn) = k2
         530    17864169 :                   exit
         531             :                end if
         532             :             end do
         533             : !
         534             : ! Identify columns with clouds in the given region.
         535             : !
         536   217832889 :             do k = k1, k2
         537   217832889 :                if (cldp(ilon,k) >= cldmin) then
         538    17864169 :                   n = n+1
         539    17864169 :                   indxmx(n,irgn) = ilon
         540    17864169 :                   exit
         541             :                endif
         542             :             end do
         543             :          endif
         544    17864169 :          ncolmx(irgn) = n
         545             : !
         546             : ! Dummy value for handling clear-sky regions
         547             : !
         548    17864169 :          indxmx(ncolmx(irgn)+1,irgn) = ncol+1
         549             : !
         550             : ! Outer loop over columns with clouds in the max-overlap region
         551             : !
         552    47172912 :          do iimx = 1, ncolmx(irgn)
         553    17864169 :             i = indxmx(iimx,irgn)
         554             : !
         555             : ! Sort cloud areas and corresponding level indices.
         556             : !
         557    17864169 :             n = 0
         558   315423093 :             do k = kx1(i,irgn),kx2(i,irgn)
         559   315423093 :                if (cldp(i,k) >= cldmin) then
         560    78579288 :                   n = n+1
         561    78579288 :                   ksort(n) = k
         562             : !
         563             : ! We need indices for clouds in order of largest to smallest, so
         564             : !    sort 1-cld in ascending order
         565             : !
         566    78579288 :                   asort(n) = 1.0_r8-cldp(i,k)
         567             :                end if
         568             :             end do
         569    17864169 :             nxs(i,irgn) = n
         570             : !
         571             : ! If nxs(i,irgn) eq 1, no need to sort.
         572             : ! If nxs(i,irgn) eq 2, sort by swapping if necessary
         573             : ! If nxs(i,irgn) ge 3, sort using local sort routine
         574             : !
         575    17864169 :             if (nxs(i,irgn) == 2) then
         576     4105690 :                if (asort(2) < asort(1)) then
         577     3130963 :                   ktmp = ksort(1)
         578     3130963 :                   ksort(1) = ksort(2)
         579     3130963 :                   ksort(2) = ktmp
         580             :
         581     3130963 :                   atmp = asort(1)
         582     3130963 :                   asort(1) = asort(2)
         583     3130963 :                   asort(2) = atmp
         584             :                endif
         585    13758479 :             else if (nxs(i,irgn) >= 3) then
         586    11268485 :                call quick_sort(asort(1:nxs(i,irgn)),ksort(1:nxs(i,irgn)))
         587             :             endif
         588             :
         589   114307626 :             do l = 1, nxs(i,irgn)
         590    96443457 :                kxs(l,i,irgn) = ksort(l)
         591             :             end do
         592             : !
         593             : ! End loop over longitude i for fluxes
         594             : !
         595             :          end do
         596             : !
         597             : ! End loop over regions irgn for max-overlap
         598             : !
         599             :       end do
         600             : !
         601             : !----------------------------------------------------------------------
         602             : ! DOWNWARD FLUXES:
         603             : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
         604             : !
         605    29308743 :       do irgn = 1, nmxrgn(ilon)
         606             : !
         607             : ! Compute clear-sky fluxes for regions without clouds
         608             : !
         609    17864169 :          iimx = 1
         610    17864169 :          if (ilon < indxmx(iimx,irgn) .and. irgn <= nrgn(ilon)) then
         611             : !
         612             : ! Calculate emissivity so that downward flux at upper boundary of region
         613             : !    can be cast in form of solution for downward flux from cloud above
         614             : !    that boundary.  Then solutions for fluxes at other levels take form of
         615             : !    random overlap expressions.  Try to locate "cloud" as close as possible
         616             : !    to TOA such that the "cloud" pseudo-emissivity is between 0 and 1.
         617             : !
         618           0 :             k1 = kx1(ilon,irgn)
         619           0 :             do km1 = ntoplw-2, k1-2
         620           0 :                km4 = km1+3
         621           0 :                k2 = k1
         622           0 :                k3 = k2+1
         623           0 :                tmp(ilon) = s(ilon,k2,min(k3,pverp))*min(1,pverp2-k3)
         624             :                emx0 = (fdl(ilon,k1)-fsdl(ilon,k1))/ &
         625           0 :                       ((fclb4(ilon,km1)-s(ilon,k2,km4)+tmp(ilon))- fsdl(ilon,k1))
         626           0 :                if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
         627             :             end do
         628           0 :             km1 = min(km1,k1-2)
         629           0 :             do k2 = kx1(ilon,irgn)+1, kx2(ilon,irgn)+1
         630           0 :                k3 = k2+1
         631           0 :                tmp(ilon) = s(ilon,k2,min(k3,pverp))*min(1,pverp2-k3)
         632             :                fdl(ilon,k2) = (1.0_r8-emx0)*fsdl(ilon,k2) + &
         633           0 :                                emx0*(fclb4(ilon,km1)-s(ilon,k2,km4)+tmp(ilon))
         634             :             end do
         635             :          else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
         636    17864169 :             iimx = iimx+1
         637             :          end if
         638             : !
         639             : ! Outer loop over columns with clouds in the max-overlap region
         640             : !
         641    47172912 :          do iimx = 1, ncolmx(irgn)
         642    17864169 :             i = indxmx(iimx,irgn)
         643             :
         644             : !
         645             : ! Calculate emissivity so that downward flux at upper boundary of region
         646             : !    can be cast in form of solution for downward flux from cloud above that
         647             : !    boundary.  Then solutions for fluxes at other levels take form of
         648             : !    random overlap expressions.  Try to locate "cloud" as close as possible
         649             : !    to TOA such that the "cloud" pseudo-emissivity is between 0 and 1.
         650             : !
         651    17864169 :             k1 = kx1(i,irgn)
         652    26075955 :             do km1 = ntoplw-2,k1-2
         653    26074827 :                km4 = km1+3
         654    26074827 :                k2 = k1
         655    26074827 :                k3 = k2 + 1
         656    26074827 :                tmp(i) = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
         657    26074827 :                tmp2(i) = s(i,k2,min(km4,pverp))*min(1,pverp2-km4)
         658    26074827 :                emx0 = (fdl(i,k1)-fsdl(i,k1))/((fclb4(i,km1)-tmp2(i)+tmp(i))-fsdl(i,k1))
         659    26075955 :                if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
         660             :             end do
         661    17864169 :             km1 = min(km1,k1-2)
         662    17864169 :             ksort(0) = km1 + 1
         663             : !
         664             : ! Loop to calculate fluxes at level k
         665             : !
         666    17864169 :             nxsk = 0
         667   333287262 :             do k = kx1(i,irgn), kx2(i,irgn)
         668             : !
         669             : ! Identify clouds (largest to smallest area) between kx1 and k
         670             : !    Since nxsk will increase with increasing k up to nxs(i,irgn), once
         671             : !    nxsk == nxs(i,irgn) then use the list constructed for previous k
         672             : !
         673   297558924 :                if (nxsk < nxs(i,irgn)) then
         674             :                   nxsk = 0
         675  1710707644 :                   do l = 1, nxs(i,irgn)
         676  1432159636 :                      k1 = kxs(l,i,irgn)
         677  1710707644 :                      if (k >= k1) then
         678   304166143 :                         nxsk = nxsk + 1
         679   304166143 :                         ksort(nxsk) = k1
         680             :                      endif
         681             :                   end do
         682             :                endif
         683             : !
         684             : ! Dummy value of index to insure computation of cloud amt is valid for l=nxsk+1
         685             : !
         686   297558924 :                ksort(nxsk+1) = pverp
         687             : !
         688             : ! Initialize iterated emissivity factors
         689             : !
         690   674792305 :                do l = 1, nxsk
         691   674792305 :                   emx(l) = emis(i,ksort(l))
         692             :                end do
         693             : !
         694             : ! Initialize iterated emissivity factor for bnd. condition at upper interface
         695             : !
         696   297558924 :                emx(0) = emx0
         697             : !
         698             : ! Initialize previous cloud amounts
         699             : !
         700   297558924 :                cld0 = 1.0_r8
         701             : !
         702             : ! Indices for flux calculations
         703             : !
         704   297558924 :                k2 = k+1
         705   297558924 :                k3 = k2+1
         706   297558924 :                tmp4 = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
         707             : !
         708             : ! Special case nxsk == 0
         709             : !
         710   297558924 :                if ( nxsk == 0 ) then
         711   199968720 :                   fdl(i,k2) = fdl(i,k2)+fsdl(i,k2)
         712   199968720 :                   if ( emx0 /= 0.0_r8 ) then
         713    29492369 :                      km1 = ksort(0)-1
         714    29492369 :                      km4 = km1+3
         715             :                      fdl(i,k2) = fdl(i,k2)+emx0* &
         716    29492369 :                                  (fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2))
         717             :                   end if
         718             :                   cycle
         719             :                end if                   ! nxsk == 0
         720             :
         721             : !
         722             : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
         723             : !
         724   572413789 :                 do l1 = 0, nxsk
         725   474823585 :                    km1 = ksort(l1)-1
         726   474823585 :                    km4 = km1+3
         727   572413789 :                    tmp3(l1) = fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2)
         728             :                end do
         729             :
         730             :                tfdl = 0.0_r8
         731             :
         732   572413789 :                do l = 1, nxsk+1
         733             : !
         734             : ! Calculate downward fluxes
         735             : !
         736   474823585 :                   cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
         737   474823585 :                   if (cld0 /= cld1) then
         738   459449201 :                      tfdl = tfdl+(cld0-cld1)*fsdl(i,k2)
         739  2249198454 :                      do l1 = 0, l - 1
         740  2249198454 :                         tfdl = tfdl+(cld0-cld1)*emx(l1)*tmp3(l1)
         741             :                      end do
         742             :                   endif
         743   474823585 :                   cld0 = cld1
         744             : !
         745             : ! Multiply emissivity factors by current cloud transmissivity
         746             : !
         747   572413789 :                   if (l <= nxsk) then
         748   377233381 :                      k1 = ksort(l)
         749   377233381 :                      trans = 1.0_r8-emis(i,k1)
         750             : !
         751             : ! Ideally the upper bound on l1 would be l-1, but the sort routine
         752             : !    scrambles the order of layers with identical cloud amounts
         753             : !
         754  3052200275 :                      do l1 = 0, nxsk
         755  3052200275 :                         if (ksort(l1) < k1) then
         756  1337483447 :                            emx(l1) = emx(l1)*trans
         757             :                         endif
         758             :                      end do
         759             :                   end if
         760             : !
         761             : ! End loop over number l of cloud levels
         762             : !
         763             :                end do
         764   115454373 :                fdl(i,k2) = tfdl
         765             : !
         766             : ! End loop over level k for fluxes
         767             : !
         768             :             end do
         769             : !
         770             : ! End loop over longitude i for fluxes
         771             : !
         772             :          end do
         773             : !
         774             : ! End loop over regions irgn for max-overlap
         775             : !
         776             :       end do
         777             :
         778             : !
         779             : !----------------------------------------------------------------------
         780             : ! UPWARD FLUXES:
         781             : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
         782             : !
         783    30057975 :       do irgn = nmxrgn(ilon), 1, -1
         784             : !
         785             : ! Compute clear-sky fluxes for regions without clouds
         786             : !
         787    17864169 :          iimx = 1
         788    17864169 :          if (ilon < indxmx(iimx,irgn) .and. irgn <= nrgn(ilon)) then
         789             : !
         790             : ! Calculate emissivity so that upward flux at lower boundary of region
         791             : !    can be cast in form of solution for upward flux from cloud below that
         792             : !    boundary.  Then solutions for fluxes at other levels take form of
         793             : !    random overlap expressions.  Try to locate "cloud" as close as possible
         794             : !    to surface such that the "cloud" pseudo-emissivity is between 0 and 1.
         795             : ! Include allowance for surface emissivity (both numerator and denominator
         796             : !    equal 1)
         797             : !
         798           0 :             k1 = kx2(ilon,irgn)+1
         799           0 :             if (k1 < pverp) then
         800           0 :                do km1 = pver-1,kx2(ilon,irgn),-1
         801           0 :                   km3 = km1+2
         802           0 :                   k2 = k1
         803           0 :                   k3 = k2+1
         804           0 :                   tmp(ilon) = s(ilon,k2,min(km3,pverp))* min(1,pverp2-km3)
         805             :                   emx0 = (ful(ilon,k1)-fsul(ilon,k1))/ &
         806           0 :                          ((fclt4(ilon,km1)+s(ilon,k2,k3)-tmp(ilon))- fsul(ilon,k1))
         807           0 :                   if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
         808             :                end do
         809           0 :                km1 = max(km1,kx2(ilon,irgn))
         810             :             else
         811           0 :                km1 = k1-1
         812           0 :                km3 = km1+2
         813           0 :                emx0 = 1.0_r8
         814             :             endif
         815             :
         816           0 :             do k2 = kx1(ilon,irgn), kx2(ilon,irgn)
         817           0 :                k3 = k2+1
         818             : !
         819             : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
         820             : !
         821           0 :                tmp(ilon) = s(ilon,k2,min(km3,pverp))* min(1,pverp2-km3)
         822             :                ful(ilon,k2) =(1.0_r8-emx0)*fsul(ilon,k2) + emx0* &
         823           0 :                              (fclt4(ilon,km1)+s(ilon,k2,k3)-tmp(ilon))
         824             :             end do
         825    17864169 :          else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
         826    17864169 :             iimx = iimx+1
         827             :          end if
         828             : !
         829             : ! Outer loop over columns with clouds in the max-overlap region
         830             : !
         831    47172912 :          do iimx = 1, ncolmx(irgn)
         832    17864169 :             i = indxmx(iimx,irgn)
         833             :
         834             : !
         835             : ! Calculate emissivity so that upward flux at lower boundary of region
         836             : !    can be cast in form of solution for upward flux from cloud at that
         837             : !    boundary.  Then solutions for fluxes at other levels take form of
         838             : !    random overlap expressions.  Try to locate "cloud" as close as possible
         839             : !    to surface such that the "cloud" pseudo-emissivity is between 0 and 1.
         840             : ! Include allowance for surface emissivity (both numerator and denominator
         841             : !    equal 1)
         842             : !
         843    17864169 :             k1 = kx2(i,irgn)+1
         844    17864169 :             if (k1 < pverp) then
         845    19704636 :                do km1 = pver-1,kx2(i,irgn),-1
         846    19702327 :                   km3 = km1+2
         847    19702327 :                   k2 = k1
         848    19702327 :                   k3 = k2+1
         849    19702327 :                   tmp(i) = s(i,k2,min(km3,pverp))*min(1,pverp2-km3)
         850    19702327 :                   emx0 = (ful(i,k1)-fsul(i,k1))/((fclt4(i,km1)+s(i,k2,k3)-tmp(i))-fsul(i,k1))
         851    19704636 :                   if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
         852             :                end do
         853     6419595 :                km1 = max(km1,kx2(i,irgn))
         854             :             else
         855             :                emx0 = 1.0_r8
         856             :                km1 = k1-1
         857             :             endif
         858    17864169 :             ksort(0) = km1 + 1
         859             :
         860             : !
         861             : ! Loop to calculate fluxes at level k
         862             : !
         863    17864169 :             nxsk = 0
         864   333287262 :             do k = kx2(i,irgn), kx1(i,irgn), -1
         865             : !
         866             : ! Identify clouds (largest to smallest area) between k and kx2
         867             : !    Since nxsk will increase with decreasing k up to nxs(i,irgn), once
         868             : !    nxsk == nxs(i,irgn) then use the list constructed for previous k
         869             : !
         870   297558924 :                if (nxsk < nxs(i,irgn)) then
         871             :                   nxsk = 0
         872   700410440 :                   do l = 1, nxs(i,irgn)
         873   602820236 :                      k1 = kxs(l,i,irgn)
         874   700410440 :                      if (k <= k1) then
         875   304166143 :                         nxsk = nxsk + 1
         876   304166143 :                         ksort(nxsk) = k1
         877             :                      endif
         878             :                   end do
         879             :                endif
         880             : !
         881             : ! Dummy value of index to insure computation of cloud amt is valid for l=nxsk+1
         882             : !
         883   297558924 :                ksort(nxsk+1) = pverp
         884             : !
         885             : ! Initialize iterated emissivity factors
         886             : !
         887  1504131705 :                do l = 1, nxsk
         888  1504131705 :                   emx(l) = emis(i,ksort(l))
         889             :                end do
         890             : !
         891             : ! Initialize iterated emissivity factor for bnd. condition at lower interface
         892             : !
         893   297558924 :                emx(0) = emx0
         894             : !
         895             : ! Initialize previous cloud amounts
         896             : !
         897   297558924 :                cld0 = 1.0_r8
         898             : !
         899             : ! Indices for flux calculations
         900             : !
         901   297558924 :                k2 = k
         902   297558924 :                k3 = k2+1
         903             : !
         904             : ! Special case nxsk == 0
         905             : !
         906   297558924 :                if ( nxsk == 0 ) then
         907    19010916 :                   ful(i,k2) = ful(i,k2)+fsul(i,k2)
         908    19010916 :                   if ( emx0 /= 0.0_r8 ) then
         909    19010916 :                      km1 = ksort(0)-1
         910    19010916 :                      km3 = km1+2
         911             : !
         912             : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
         913             : !
         914             :                      ful(i,k2) = ful(i,k2)+emx0* &
         915    19010916 :                                 (fclt4(i,km1)+s(i,k2,k3)-(s(i,k2,min(km3,pverp))*min(1,pverp2-km3))-fsul(i,k2))
         916             :
         917             :                   end if
         918             :                   cycle
         919             :                end if
         920             : !
         921             : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
         922             : !
         923  1763668797 :                do l1 = 0, nxsk
         924  1485120789 :                   km1 = ksort(l1)-1
         925  1485120789 :                   km3 = km1+2
         926             : !
         927             : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
         928             : !
         929  1763668797 :                   tmp3(l1) = fclt4(i,km1)+s(i,k2,k3)-(s(i,k2,min(km3,pverp))*min(1,pverp2-km3))- fsul(i,k2)
         930             :                end do
         931             :
         932             :                tful = 0.0_r8
         933             :
         934  1763668797 :                do l = 1, nxsk+1
         935             : !
         936             : ! Calculate upward fluxes
         937             : !
         938  1485120789 :                   cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
         939  1485120789 :                   if (cld0 /= cld1) then
         940  1449201158 :                      tful = tful+(cld0-cld1)*fsul(i,k2)
         941  7420388250 :                      do l1 = 0, l - 1
         942  7420388250 :                         tful = tful+(cld0-cld1)*emx(l1)*tmp3(l1)
         943             :                      end do
         944             :                   endif
         945  1485120789 :                   cld0 = cld1
         946             : !
         947             : ! Multiply emissivity factors by current cloud transmissivity
         948             : !
         949  1763668797 :                   if (l <= nxsk) then
         950  1206572781 :                      k1 = ksort(l)
         951  1206572781 :                      trans = 1.0_r8-emis(i,k1)
         952             : !
         953             : ! Ideally the upper bound on l1 would be l-1, but the sort routine
         954             : !    scrambles the order of layers with identical cloud amounts
         955             : !
         956 10275808737 :                      do l1 = 0, nxsk
         957 10275808737 :                         if (ksort(l1) > k1) then
         958  4534617978 :                            emx(l1) = emx(l1)*trans
         959             :                         endif
         960             :                      end do
         961             :                   end if
         962             : !
         963             : ! End loop over number l of cloud levels
         964             : !
         965             :                end do
         966   296412177 :                ful(i,k2) = tful
         967             : !
         968             : ! End loop over level k for fluxes
         969             : !
         970             :             end do
         971             : !
         972             : ! End loop over longitude i for fluxes
         973             : !
         974             :          end do
         975             : !
         976             : ! End loop over regions irgn for max-overlap
         977             : !
         978             :       end do
         979             : !
         980             : ! End outermost longitude loop
         981             : !
         982             :    end do
         983             : !CSD$ END PARALLEL DO
         984             : !
         985             : ! End cloud modification loops
         986             : !
         987             : !----------------------------------------------------------------------
         988             : ! All longitudes: store history tape quantities
         989             : !
         990    12510432 :    do i=1,ncol
         991    11761200 :       flwds(i) = fdl (i,pverp )
         992    11761200 :       fldsc(i) = fsdl(i,pverp )
         993    11761200 :       flns(i)  = ful (i,pverp ) - fdl (i,pverp )
         994    11761200 :       flnsc(i) = fsul(i,pverp ) - fsdl(i,pverp )
         995    11761200 :       flnt(i)  = ful (i,ntoplw) - fdl (i,ntoplw)
         996    11761200 :       flntc(i) = fsul(i,ntoplw) - fsdl(i,ntoplw)
         997    11761200 :       flut(i)  = ful (i,ntoplw)
         998    12510432 :       flutc(i) = fsul(i,ntoplw)
         999             :    end do
        1000             :
        1001      749232 :    if (single_column.and.scm_crm_mode) then
        1002           0 :       call outfld('FUL     ',ful*1.e-3_r8,pcols,lchnk)
        1003           0 :       call outfld('FDL     ',fdl*1.e-3_r8,pcols,lchnk)
        1004           0 :       call outfld('FULC    ',fsul*1.e-3_r8,pcols,lchnk)
        1005           0 :       call outfld('FDLC    ',fsdl*1.e-3_r8,pcols,lchnk)
        1006             :    endif
        1007             :
        1008    20978496 :    do k = ntoplw,pverp
        1009   338530896 :       do i = 1,ncol
        1010   337781664 :          fnl(i,k) = ful(i,k) - fdl(i,k)
        1011             :       end do
        1012             :    end do
        1013             : !
        1014             : ! Computation of longwave heating (J/kg/s)
        1015             : !
        1016    20229264 :    do k=ntoplw,pver
        1017   326020464 :       do i=1,ncol
        1018   917373600 :          qrl(i,k) = (ful(i,k) - fdl(i,k) - ful(i,k+1) + fdl(i,k+1))* &
        1019   917373600 :                      1.E-4_r8*gravit_cgs/((pint(i,k) - pint(i,k+1)))
        1020             :          qrlc(i,k) = (fsul(i,k) - fsdl(i,k) - fsul(i,k+1) + fsdl(i,k+1))* &
        1021   325271232 :                       1.E-4_r8*gravit_cgs/((pint(i,k) - pint(i,k+1)))
        1022             :       end do
        1023             :    end do
        1024             : ! Return 0 above solution domain
        1025      749232 :    if ( ntoplw > 1 )then
        1026           0 :       qrl(:ncol,:ntoplw-1) = 0._r8
        1027           0 :       qrlc(:ncol,:ntoplw-1) = 0._r8
        1028             :    end if
        1029             : !
        1030      749232 :    return
        1031     1498464 : end subroutine radclwmx
        1032             :
        1033             :
        1034             :
        1035             : !-------------------------------------------------------------------------------
        1036             :
        1037        1536 : subroutine radlw_init(gravit, stebol)
        1038             : !-----------------------------------------------------------------------
        1039             : !
        1040             : ! Purpose:
        1041             : ! Initialize various constants for radiation scheme; note that
        1042             : ! the radiation scheme uses cgs units.
        1043             : !-----------------------------------------------------------------------
        1044             : !
        1045             : ! Input arguments
        1046             : !
        1047             :    real(r8), intent(in) :: gravit     ! Acceleration of gravity (MKS)
        1048             :    real(r8), intent(in) :: stebol     ! Stefan-Boltzmann's constant (MKS)
        1049             : !
        1050             : !-----------------------------------------------------------------------
        1051             : !
        1052             : ! Set general radiation consts; convert to cgs units where appropriate:
        1053             : !
        1054        1536 :    gravit_cgs  =  100._r8*gravit
        1055        1536 :    stebol_cgs  =  1.e3_r8*stebol
        1056             :
        1057      749232 : end subroutine radlw_init
        1058             :
        1059             : !-------------------------------------------------------------------------------
        1060       68112 : subroutine aer_trans_from_od(ncol, odap_aer, aer_trn_ttl)
        1061             :   use radconstants,  only: nlwbands
        1062             :   use ppgrid,        only: pcols, pver, pverp
        1063             :   integer, intent(in) :: ncol
        1064             : ! [fraction] absorption optical depth, per layer
        1065             :   real(r8), intent(in) :: odap_aer(pcols,pver,nlwbands)
        1066             : ! [fraction] Total transmission between interfaces k1 and k2
        1067             :   real(r8), intent(out) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands)
        1068             : ! [fraction] absorption optical depth, cumulative from top
        1069             :   real(r8) :: odap_aer_ttl(pcols,pverp,nlwbands)
        1070             :
        1071             :   integer i, k1, k2, bnd_idx ! column iterator, level iterators, band iterator
        1072             :   ! odap_aer_ttl is cumulative total optical depth from top of atmosphere
        1073       68112 :   odap_aer_ttl = 0._r8
        1074      544896 :   do bnd_idx=1,nlwbands
        1075    12941280 :     do k1=1,pver
        1076   207467568 :       do i=1,ncol
        1077   206990784 :         odap_aer_ttl(i,k1+1,bnd_idx) = odap_aer_ttl(i,k1,bnd_idx) +  odap_aer(i,k1,bnd_idx)
        1078             :       end do
        1079             :     end do
        1080             :   end do
        1081             :
        1082             :   ! compute transmission from top of atmosphere to level
        1083             :   ! where angular dependence of optical depth has been
        1084             :   ! integrated (giving factor 1.666666)
        1085     1907136 :   do k1=1,pverp
        1086   220750992 :      aer_trn_ttl(1:pcols,k1,k1,:)=1._r8
        1087             :   enddo
        1088             :   aer_trn_ttl(1:ncol,1,2:pverp,1:nlwbands) = &
        1089   207535680 :     exp(-1.66_r8 * odap_aer_ttl(1:ncol,2:pverp,1:nlwbands) )
        1090             :
        1091             :   ! compute transmission between a given layer (k1) and lower layers.
        1092     1770912 :   do k1=2,pver
        1093    23907312 :     do k2=k1+1,pverp
        1094             :       aer_trn_ttl(1:ncol,k1,k2,1:nlwbands) = &
        1095    22136400 :           aer_trn_ttl(1:ncol,1,k2,1:nlwbands) / &
        1096  2633360400 :           aer_trn_ttl(1:ncol,1,k1,1:nlwbands)
        1097             :     end do
        1098             :   end do
        1099             :
        1100             :   ! transmission from k1 to k2 is same as transmission from k2 to k1.
        1101     1839024 :   do k1=2,pverp
        1102    25746336 :     do k2=1,k1-1
        1103  2820053808 :         aer_trn_ttl(1:ncol,k1,k2,1:nlwbands)=aer_trn_ttl(1:ncol,k2,k1,1:nlwbands)
        1104             :     end do
        1105             :   end do
        1106             :
        1107       68112 :   return
        1108             : end subroutine aer_trans_from_od
        1109             :
        1110             : end module radlw

Generated by: LCOV version 1.14