LCOV - code coverage report
Current view: | top level - physics/rrtmgp - radiation.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 950 | 1061 | 89.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 23 | 24 | 95.8 %

              Line data    Source code
           1             : module radiation
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : !
           5             : ! CAM interface to RRTMGP radiation parameterization.
           6             : !
           7             : !---------------------------------------------------------------------------------
           8             :
           9             : use shr_kind_mod,        only: r8=>shr_kind_r8, cl=>shr_kind_cl
          10             : use spmd_utils,          only: masterproc
          11             : use ppgrid,              only: pcols, pver, pverp, begchunk, endchunk
          12             : use ref_pres,            only: pref_edge
          13             : use physics_types,       only: physics_state, physics_ptend
          14             : use phys_control,        only: phys_getopts
          15             : use physics_buffer,      only: physics_buffer_desc, pbuf_add_field, dtype_r8, pbuf_get_index, &
          16             :                                pbuf_set_field, pbuf_get_field, pbuf_old_tim_idx
          17             : use camsrfexch,          only: cam_out_t, cam_in_t
          18             : use physconst,           only: cappa, cpair, gravit
          19             : use solar_irrad_data,    only: sol_tsi
          20             :
          21             : use time_manager,        only: get_nstep, is_first_step, is_first_restart_step, &
          22             :                                get_curr_calday, get_step_size
          23             :
          24             : use rad_constituents,    only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_gas, rad_cnst_out
          25             :
          26             : use rrtmgp_inputs,       only: rrtmgp_inputs_init
          27             :
          28             : use radconstants,        only: nradgas, gasnamelength, gaslist, nswbands, nlwbands, &
          29             :                                nswgpts, set_wavenumber_bands
          30             :
          31             : use cloud_rad_props,     only: cloud_rad_props_init
          32             :
          33             : use cospsimulator_intr,  only: docosp, cospsimulator_intr_init, &
          34             :                                cospsimulator_intr_run, cosp_nradsteps
          35             :
          36             : use scamMod,             only: scm_crm_mode, single_column, have_cld, cldobs
          37             :
          38             : use cam_history,         only: addfld, add_default, horiz_only, outfld, hist_fld_active
          39             :
          40             : use radiation_data,      only: rad_data_register, rad_data_init
          41             :
          42             : use ioFileMod,           only: getfil
          43             : use cam_pio_utils,       only: cam_pio_openfile
          44             : use pio,                 only: file_desc_t, var_desc_t,                       &
          45             :                                pio_int, pio_double, PIO_NOERR,                &
          46             :                                pio_seterrorhandling, PIO_BCAST_ERROR,         &
          47             :                                pio_inq_dimlen, pio_inq_dimid, pio_inq_varid,  &
          48             :                                pio_def_var, pio_put_var, pio_get_var,         &
          49             :                                pio_put_att, PIO_NOWRITE, pio_closefile
          50             :
          51             : use mo_gas_concentrations, only: ty_gas_concs
          52             : use mo_gas_optics_rrtmgp,  only: ty_gas_optics_rrtmgp
          53             : use mo_optical_props,      only: ty_optical_props_1scl, ty_optical_props_2str
          54             : use mo_source_functions,   only: ty_source_func_lw
          55             : use mo_fluxes,             only: ty_fluxes_broadband
          56             : use mo_fluxes_byband,      only: ty_fluxes_byband
          57             :
          58             : use string_utils,        only: to_lower
          59             : use cam_abortutils,      only: endrun, handle_allocate_error
          60             : use cam_logfile,         only: iulog
          61             :
          62             :
          63             : implicit none
          64             : private
          65             : save
          66             :
          67             : public :: &
          68             :    radiation_readnl,         &! read namelist variables
          69             :    radiation_register,       &! registers radiation physics buffer fields
          70             :    radiation_do,             &! query which radiation calcs are done this timestep
          71             :    radiation_init,           &! initialization
          72             :    radiation_define_restart, &! define variables for restart
          73             :    radiation_write_restart,  &! write variables to restart
          74             :    radiation_read_restart,   &! read variables from restart
          75             :    radiation_tend,           &! compute heating rates and fluxes
          76             :    rad_out_t                  ! type for diagnostic outputs
          77             :
          78             : integer,public, allocatable :: cosp_cnt(:)       ! counter for cosp
          79             : integer,public              :: cosp_cnt_init = 0 !initial value for cosp counter
          80             :
          81             : real(r8), public, protected :: nextsw_cday       ! future radiation calday for surface models
          82             :
          83             : type rad_out_t
          84             :    real(r8) :: solin(pcols)         ! Solar incident flux
          85             :
          86             :    real(r8) :: qrsc(pcols,pver)
          87             :
          88             :    real(r8) :: fsnsc(pcols)         ! Clear sky surface abs solar flux
          89             :    real(r8) :: fsntc(pcols)         ! Clear sky total column abs solar flux
          90             :    real(r8) :: fsdsc(pcols)         ! Clear sky surface downwelling solar flux
          91             :
          92             :    real(r8) :: fsntoa(pcols)        ! Net solar flux at TOA
          93             :    real(r8) :: fsntoac(pcols)       ! Clear sky net solar flux at TOA
          94             :    real(r8) :: fsutoa(pcols)        ! upwelling solar flux at TOA
          95             :
          96             :    real(r8) :: fsnirt(pcols)        ! Near-IR flux absorbed at toa
          97             :    real(r8) :: fsnrtc(pcols)        ! Clear sky near-IR flux absorbed at toa
          98             :    real(r8) :: fsnirtsq(pcols)      ! Near-IR flux absorbed at toa >= 0.7 microns
          99             :
         100             :    real(r8) :: fsn200(pcols)        ! Net SW flux interpolated to 200 mb
         101             :    real(r8) :: fsn200c(pcols)       ! Net clear-sky SW flux interpolated to 200 mb
         102             :    real(r8) :: fsnr(pcols)          ! Net SW flux interpolated to tropopause
         103             :
         104             :    real(r8) :: flux_sw_up(pcols,pverp)     ! upward shortwave flux on interfaces
         105             :    real(r8) :: flux_sw_clr_up(pcols,pverp) ! upward shortwave clearsky flux
         106             :    real(r8) :: flux_sw_dn(pcols,pverp)     ! downward flux
         107             :    real(r8) :: flux_sw_clr_dn(pcols,pverp) ! downward clearsky flux
         108             :
         109             :    real(r8) :: flux_lw_up(pcols,pverp)     ! upward longwave flux on interfaces
         110             :    real(r8) :: flux_lw_clr_up(pcols,pverp) ! upward longwave clearsky flux
         111             :    real(r8) :: flux_lw_dn(pcols,pverp)     ! downward flux
         112             :    real(r8) :: flux_lw_clr_dn(pcols,pverp) ! downward clearsky flux
         113             :
         114             :    real(r8) :: qrlc(pcols,pver)
         115             :
         116             :    real(r8) :: flntc(pcols)         ! Clear sky lw flux at model top
         117             :    real(r8) :: flut(pcols)          ! Upward flux at top of model
         118             :    real(r8) :: flutc(pcols)         ! Upward Clear Sky flux at top of model
         119             :    real(r8) :: lwcf(pcols)          ! longwave cloud forcing
         120             :
         121             :    real(r8) :: fln200(pcols)        ! net longwave flux interpolated to 200 mb
         122             :    real(r8) :: fln200c(pcols)       ! net clearsky longwave flux interpolated to 200 mb
         123             :    real(r8) :: flnr(pcols)          ! net longwave flux interpolated to tropopause
         124             :
         125             :    real(r8) :: flnsc(pcols)         ! Clear sky lw flux at srf (up-down)
         126             :    real(r8) :: fldsc(pcols)         ! Clear sky lw flux at srf (down)
         127             :
         128             :    real(r8) :: tot_cld_vistau(pcols,pver)   ! gbx water+ice cloud optical depth (only during day, night = fillvalue)
         129             :    real(r8) :: tot_icld_vistau(pcols,pver)  ! in-cld water+ice cloud optical depth (only during day, night = fillvalue)
         130             :    real(r8) :: liq_icld_vistau(pcols,pver)  ! in-cld liq cloud optical depth (only during day, night = fillvalue)
         131             :    real(r8) :: ice_icld_vistau(pcols,pver)  ! in-cld ice cloud optical depth (only during day, night = fillvalue)
         132             :    real(r8) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth for output on history files
         133             :    real(r8) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth for output on history files
         134             : end type rad_out_t
         135             :
         136             : ! Control variables set via namelist
         137             : character(len=cl) :: coefs_lw_file ! filepath for lw coefficients
         138             : character(len=cl) :: coefs_sw_file ! filepath for sw coefficients
         139             :
         140             : integer :: iradsw = -1     ! freq. of shortwave radiation calc in time steps (positive)
         141             :                            ! or hours (negative).
         142             : integer :: iradlw = -1     ! frequency of longwave rad. calc. in time steps (positive)
         143             :                            ! or hours (negative).
         144             :
         145             : integer :: irad_always = 0 ! Specifies length of time in timesteps (positive)
         146             :                            ! or hours (negative) SW/LW radiation will be
         147             :                            ! run continuously from the start of an
         148             :                            ! initial or restart run
         149             : logical :: use_rad_dt_cosz  = .false. ! if true, use radiation dt for all cosz calculations
         150             : logical :: spectralflux     = .false. ! calculate fluxes (up and down) per band.
         151             : logical :: graupel_in_rad   = .false. ! graupel in radiation code
         152             : logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation
         153             :
         154             : ! active_calls is set by a rad_constituents method after parsing namelist input
         155             : ! for the rad_climate and rad_diag_N entries.
         156             : logical :: active_calls(0:N_DIAG)
         157             :
         158             : ! Physics buffer indices
         159             : integer :: qrs_idx      = 0
         160             : integer :: qrl_idx      = 0
         161             : integer :: su_idx       = 0
         162             : integer :: sd_idx       = 0
         163             : integer :: lu_idx       = 0
         164             : integer :: ld_idx       = 0
         165             : integer :: fsds_idx     = 0
         166             : integer :: fsns_idx     = 0
         167             : integer :: fsnt_idx     = 0
         168             : integer :: flns_idx     = 0
         169             : integer :: flnt_idx     = 0
         170             : integer :: cld_idx      = 0
         171             : integer :: cldfsnow_idx = 0
         172             : integer :: cldfgrau_idx = 0
         173             :
         174             : character(len=4) :: diag(0:N_DIAG) =(/'    ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ',&
         175             :                                       '_d6 ','_d7 ','_d8 ','_d9 ','_d10'/)
         176             :
         177             : ! averaging time interval for zenith angle
         178             : real(r8) :: dt_avg = 0._r8
         179             : real(r8) :: rad_uniform_angle = -99._r8
         180             :
         181             : ! Number of layers in radiation calculations.
         182             : integer :: nlay
         183             :
         184             : ! Number of CAM layers in radiation calculations.  Is either equal to nlay, or is
         185             : ! 1 less than nlay if "extra layer" is used in the radiation calculations.
         186             : integer :: nlaycam
         187             :
         188             : ! Indices for copying data between CAM/WACCM and RRTMGP arrays.  Since RRTMGP is
         189             : ! vertical order agnostic we can send data using the top to bottom order used
         190             : ! in CAM/WACCM.  But the number of layers that RRTMGP does computations for
         191             : ! may not match the number of layers in CAM/WACCM for two reasons:
         192             : ! 1. If the CAM model top is below 1 Pa, then RRTMGP does calculations for an
         193             : !    extra layer that is added between 1 Pa and the model top.
         194             : ! 2. If the WACCM model top is above 1 Pa, then RRMTGP only does calculations
         195             : !    for those model layers that are below 1 Pa.
         196             : integer :: ktopcam ! Index in CAM arrays of top level (layer or interface) at which
         197             :                    ! RRTMGP is active.
         198             : integer :: ktoprad ! Index in RRTMGP arrays of the layer or interface corresponding
         199             :                    ! to CAM's top layer or interface.
         200             :                    ! Note: for CAM's top to bottom indexing, the index of a given layer
         201             :                    ! (midpoint) and the upper interface of that layer, are the same.
         202             :
         203             : ! Gas optics objects contain the data read from the coefficients files.
         204             : type(ty_gas_optics_rrtmgp) :: kdist_sw
         205             : type(ty_gas_optics_rrtmgp) :: kdist_lw
         206             :
         207             : ! lower case version of gaslist for RRTMGP
         208             : character(len=gasnamelength) :: gaslist_lc(nradgas)
         209             :
         210             : type(var_desc_t) :: cospcnt_desc  ! cosp
         211             : type(var_desc_t) :: nextsw_cday_desc
         212             :
         213             : !=========================================================================================
         214             : contains
         215             : !=========================================================================================
         216             :
         217        1536 : subroutine radiation_readnl(nlfile)
         218             :
         219             :    ! Read radiation_nl namelist group.
         220             :
         221             :    use namelist_utils,  only: find_group_name
         222             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, &
         223             :                               mpi_character, mpi_real8
         224             :
         225             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         226             :
         227             :    ! Local variables
         228             :    integer :: unitn, ierr
         229             :    integer :: dtime      ! timestep size
         230             :    character(len=32) :: errmsg
         231             :    character(len=*), parameter :: sub = 'radiation_readnl'
         232             :
         233             :    character(len=cl) :: rrtmgp_coefs_lw_file, rrtmgp_coefs_sw_file
         234             :
         235             :    namelist /radiation_nl/ &
         236             :       rrtmgp_coefs_lw_file, rrtmgp_coefs_sw_file, iradsw, iradlw,        &
         237             :       irad_always, use_rad_dt_cosz, spectralflux, use_rad_uniform_angle, &
         238             :       rad_uniform_angle, graupel_in_rad
         239             :    !-----------------------------------------------------------------------------
         240             :
         241        1536 :    if (masterproc) then
         242           2 :       open( newunit=unitn, file=trim(nlfile), status='old' )
         243           2 :       call find_group_name(unitn, 'radiation_nl', status=ierr)
         244           2 :       if (ierr == 0) then
         245           2 :          read(unitn, radiation_nl, iostat=ierr)
         246           2 :          if (ierr /= 0) then
         247           0 :             write(errmsg,'(a,i5)') 'iostat =', ierr
         248           0 :             call endrun(sub//': ERROR reading namelist: '//trim(errmsg))
         249             :          end if
         250             :       end if
         251           2 :       close(unitn)
         252             :    end if
         253             :
         254             :    ! Broadcast namelist variables
         255        1536 :    call mpi_bcast(rrtmgp_coefs_lw_file, cl, mpi_character, mstrid, mpicom, ierr)
         256        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rrtmgp_coefs_lw_file")
         257        1536 :    call mpi_bcast(rrtmgp_coefs_sw_file, cl, mpi_character, mstrid, mpicom, ierr)
         258        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rrtmgp_coefs_sw_file")
         259        1536 :    call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr)
         260        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw")
         261        1536 :    call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr)
         262        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw")
         263        1536 :    call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr)
         264        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always")
         265        1536 :    call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr)
         266        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz")
         267        1536 :    call mpi_bcast(spectralflux, 1, mpi_logical, mstrid, mpicom, ierr)
         268        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: spectralflux")
         269        1536 :    call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr)
         270        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle")
         271        1536 :    call mpi_bcast(rad_uniform_angle, 1, mpi_real8, mstrid, mpicom, ierr)
         272        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle")
         273        1536 :    call mpi_bcast(graupel_in_rad, 1, mpi_logical, mstrid, mpicom, ierr)
         274        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: graupel_in_rad")
         275             :
         276        1536 :    if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then
         277             :       call endrun(sub//': ERROR - use_rad_uniform_angle is set to .true,' &
         278           0 :                      //' but rad_uniform_angle is not set ')
         279             :    end if
         280             :
         281             :    ! Set module data
         282        1536 :    coefs_lw_file   = rrtmgp_coefs_lw_file
         283        1536 :    coefs_sw_file   = rrtmgp_coefs_sw_file
         284             :
         285             :    ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
         286        1536 :    dtime  = get_step_size()
         287        1536 :    if (iradsw      < 0) iradsw      = nint((-iradsw     *3600._r8)/dtime)
         288        1536 :    if (iradlw      < 0) iradlw      = nint((-iradlw     *3600._r8)/dtime)
         289        1536 :    if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime)
         290             :
         291             :    !-----------------------------------------------------------------------
         292             :    ! Print runtime options to log.
         293             :    !-----------------------------------------------------------------------
         294             :
         295        1536 :    if (masterproc) then
         296           2 :       write(iulog,*) 'RRTMGP radiation scheme parameters:'
         297           2 :       write(iulog,10) trim(coefs_lw_file), trim(coefs_sw_file), nlwbands, nswbands, &
         298           4 :          iradsw, iradlw, irad_always, use_rad_dt_cosz, spectralflux, graupel_in_rad
         299             :    end if
         300             :
         301             : 10 format('  LW coefficents file: ',                                a/, &
         302             :           '  SW coefficents file: ',                                a/, &
         303             :           '  Number of LW bands:                                 ',i5/, &
         304             :           '  Number of SW bands:                                 ',i5/, &
         305             :           '  Frequency (timesteps) of Shortwave Radiation calc:  ',i5/, &
         306             :           '  Frequency (timesteps) of Longwave Radiation calc:   ',i5/, &
         307             :           '  SW/LW calc done every timestep for first N steps. N=',i5/, &
         308             :           '  Use average zenith angle:                           ',l5/, &
         309             :           '  Output spectrally resolved fluxes:                  ',l5/, &
         310             :           '  Graupel in Radiation Code:                          ',l5/)
         311             :
         312        1536 : end subroutine radiation_readnl
         313             :
         314             : !================================================================================================
         315             :
         316        1536 : subroutine radiation_register
         317             :
         318             :    ! Register radiation fields in the physics buffer
         319             :
         320        1536 :    call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate
         321        1536 :    call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave  radiative heating rate
         322             :
         323        1536 :    call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux
         324        1536 :    call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux
         325        1536 :    call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux
         326             :
         327        1536 :    call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux
         328        1536 :    call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux
         329             :
         330             :    ! If the namelist has been configured for preserving the spectral fluxes, then create
         331             :    ! physics buffer variables to store the results.  This data is accessed by CARMA.
         332        1536 :    if (spectralflux) then
         333           0 :       call pbuf_add_field('SU'  , 'global',dtype_r8,(/pcols,pverp,nswbands/), su_idx) ! shortwave upward flux (per band)
         334           0 :       call pbuf_add_field('SD'  , 'global',dtype_r8,(/pcols,pverp,nswbands/), sd_idx) ! shortwave downward flux (per band)
         335           0 :       call pbuf_add_field('LU'  , 'global',dtype_r8,(/pcols,pverp,nlwbands/), lu_idx) ! longwave upward flux (per band)
         336           0 :       call pbuf_add_field('LD'  , 'global',dtype_r8,(/pcols,pverp,nlwbands/), ld_idx) ! longwave downward flux (per band)
         337             :    end if
         338             :
         339             :    ! Register fields for offline radiation driver.
         340        1536 :    call rad_data_register()
         341             :
         342        1536 : end subroutine radiation_register
         343             :
         344             : !================================================================================================
         345             :
         346     5232240 : function radiation_do(op, timestep)
         347             :
         348             :    ! Return true if the specified operation is done this timestep.
         349             :
         350             :    character(len=*), intent(in) :: op             ! name of operation
         351             :    integer, intent(in), optional:: timestep
         352             :    logical                      :: radiation_do   ! return value
         353             :
         354             :    ! Local variables
         355             :    integer :: nstep             ! current timestep number
         356             :    !-----------------------------------------------------------------------
         357             :
         358     5232240 :    if (present(timestep)) then
         359     5232240 :       nstep = timestep
         360             :    else
         361           0 :       nstep = get_nstep()
         362             :    end if
         363             :
         364     3739968 :    select case (op)
         365             :       case ('sw') ! do a shortwave heating calc this timestep?
         366             :          radiation_do = nstep == 0  .or.  iradsw == 1                     &
         367             :                      .or. (mod(nstep-1,iradsw) == 0  .and.  nstep /= 1)   &
         368     3739968 :                      .or. nstep <= irad_always
         369             :       case ('lw') ! do a longwave heating calc this timestep?
         370             :          radiation_do = nstep == 0  .or.  iradlw == 1                     &
         371             :                      .or. (mod(nstep-1,iradlw) == 0  .and.  nstep /= 1)   &
         372     1492272 :                      .or. nstep <= irad_always
         373             :       case default
         374     5232240 :          call endrun('radiation_do: unknown operation:'//op)
         375             :    end select
         376             :
         377     5232240 : end function radiation_do
         378             :
         379             : !================================================================================================
         380             :
         381     1492272 : real(r8) function radiation_nextsw_cday()
         382             :
         383             :    ! If a SW radiation calculation will be done on the next time-step, then return
         384             :    ! the calendar day of that time-step.  Otherwise return -1.0
         385             :
         386             :    ! Local variables
         387             :    integer :: nstep      ! timestep counter
         388             :    logical :: dosw       ! true => do shosrtwave calc
         389             :    integer :: offset     ! offset for calendar day calculation
         390             :    integer :: dtime      ! integer timestep size
         391             :    real(r8):: caldayp1   ! calendar day of next time-step
         392             :
         393             :    !-----------------------------------------------------------------------
         394             :
         395     1492272 :    radiation_nextsw_cday = -1._r8
         396     1492272 :    dosw   = .false.
         397     1492272 :    nstep  = get_nstep()
         398     1492272 :    dtime  = get_step_size()
         399     1492272 :    offset = 0
         400             :    do while (.not. dosw)
         401     2247696 :       nstep = nstep + 1
         402     2247696 :       offset = offset + dtime
         403     2247696 :       if (radiation_do('sw', nstep)) then
         404     1492272 :          radiation_nextsw_cday = get_curr_calday(offset=offset)
         405             :          dosw = .true.
         406             :       end if
         407             :    end do
         408     1492272 :    if(radiation_nextsw_cday == -1._r8) then
         409           0 :       call endrun('error in radiation_nextsw_cday')
         410             :    end if
         411             :
         412             :    ! determine if next radiation time-step not equal to next time-step
         413     1492272 :    if (get_nstep() >= 1) then
         414     1486080 :       caldayp1 = get_curr_calday(offset=int(dtime))
         415     1486080 :       if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
         416             :    end if
         417             :
         418     1492272 : end function radiation_nextsw_cday
         419             :
         420             : !================================================================================================
         421             :
         422        1536 : subroutine radiation_init(pbuf2d)
         423             :
         424             :    ! Initialize the radiation and cloud optics.
         425             :    ! Add fields to the history buffer.
         426             :
         427             :    ! arguments
         428             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         429             :
         430             :    ! local variables
         431             :    character(len=128) :: errmsg
         432             :
         433             :    ! names of gases that are available in the model
         434             :    ! -- needed for the kdist initialization routines
         435        1536 :    type(ty_gas_concs) :: available_gases
         436             :
         437             :    integer :: i, icall
         438             :    integer :: nstep                       ! current timestep number
         439             :    logical :: history_amwg                ! output the variables used by the AMWG diag package
         440             :    logical :: history_vdiag               ! output the variables used by the AMWG variability diag package
         441             :    logical :: history_budget              ! output tendencies and state variables for CAM4
         442             :                                           ! temperature, water vapor, cloud ice and cloud
         443             :                                           ! liquid budgets.
         444             :    integer :: history_budget_histfile_num ! history file number for budget fields
         445             :    integer :: ierr, istat
         446             :
         447             :    integer :: dtime
         448             :
         449             :    character(len=*), parameter :: sub = 'radiation_init'
         450             :    !-----------------------------------------------------------------------
         451             :
         452             :    ! Number of layers in radiation calculation is capped by the number of
         453             :    ! pressure interfaces below 1 Pa.  When the entire model atmosphere is
         454             :    ! below 1 Pa then an extra layer is added to the top of the model for
         455             :    ! the purpose of the radiation calculation.
         456             :
         457      145920 :    nlay = count( pref_edge(:) > 1._r8 ) ! pascals (0.01 mbar)
         458             :
         459        1536 :    if (nlay == pverp) then
         460             :       ! Top model interface is below 1 Pa.  RRTMGP is active in all model layers plus
         461             :       ! 1 extra layer between model top and 1 Pa.
         462           0 :       ktopcam = 1
         463           0 :       ktoprad = 2
         464           0 :       nlaycam = pver
         465        1536 :    else if (nlay == (pverp-1)) then
         466             :       ! Special case nlay == (pverp-1) -- topmost interface outside bounds (CAM MT config), treat as if it is ok.
         467        1536 :       ktopcam = 1
         468        1536 :       ktoprad = 2
         469        1536 :       nlaycam = pver
         470        1536 :       nlay = nlay+1 ! reassign the value so later code understands to treat this case like nlay==pverp
         471        1536 :       write(iulog,*) 'RADIATION_INIT: Special case of 1 model interface at p < 1Pa. Top layer will be INCLUDED in radiation calculation.'
         472        1536 :       write(iulog,*) 'RADIATION_INIT: nlay = ',nlay, ' same as pverp: ',nlay==pverp
         473             :    else
         474             :       ! nlay < pverp.  nlay layers are used in radiation calcs, and they are
         475             :       ! all CAM layers.
         476           0 :       ktopcam = pver - nlay + 1
         477           0 :       ktoprad = 1
         478           0 :       nlaycam = nlay
         479             :    end if
         480             :
         481             :    ! Create lowercase version of the gaslist for RRTMGP.  The ty_gas_concs objects
         482             :    ! work with CAM's uppercase names, but other objects that get input from the gas
         483             :    ! concs objects don't work.
         484       13824 :    do i = 1, nradgas
         485       13824 :       gaslist_lc(i) = to_lower(gaslist(i))
         486             :    end do
         487             :
         488        1536 :    errmsg = available_gases%init(gaslist_lc)
         489        1536 :    call stop_on_err(errmsg, sub, 'available_gases%init')
         490             :
         491             :    ! Read RRTMGP coefficients files and initialize kdist objects.
         492        1536 :    call coefs_init(coefs_sw_file, available_gases, kdist_sw)
         493        1536 :    call coefs_init(coefs_lw_file, available_gases, kdist_lw)
         494             :
         495             :    ! Set the sw/lw band boundaries in radconstants.  Also sets
         496             :    ! indicies of specific bands for diagnostic output and COSP input.
         497        1536 :    call set_wavenumber_bands(kdist_sw, kdist_lw)
         498             :
         499             :    ! The spectral band boundaries need to be set before this init is called.
         500        1536 :    call rrtmgp_inputs_init(ktopcam, ktoprad)
         501             :
         502             :    ! initialize output fields for offline driver
         503        1536 :    call rad_data_init(pbuf2d)
         504             :
         505        1536 :    call cloud_rad_props_init()
         506             :
         507        1536 :    cld_idx      = pbuf_get_index('CLD')
         508        1536 :    cldfsnow_idx = pbuf_get_index('CLDFSNOW', errcode=ierr)
         509        1536 :    cldfgrau_idx = pbuf_get_index('CLDFGRAU', errcode=ierr)
         510             :
         511        1536 :    if (is_first_step()) then
         512         768 :       call pbuf_set_field(pbuf2d, qrl_idx, 0._r8)
         513             :    end if
         514             :
         515             :    ! Set the radiation timestep for cosz calculations if requested using
         516             :    ! the adjusted iradsw value from radiation
         517        1536 :    if (use_rad_dt_cosz)  then
         518           0 :       dtime  = get_step_size()
         519           0 :       dt_avg = iradsw*dtime
         520             :    end if
         521             :
         522             :    ! Surface components to get radiation computed today
         523        1536 :    if (.not. is_first_restart_step()) then
         524         768 :       nextsw_cday = get_curr_calday()
         525             :    end if
         526             :
         527             :    call phys_getopts(history_amwg_out   = history_amwg,    &
         528             :                      history_vdiag_out  = history_vdiag,   &
         529             :                      history_budget_out = history_budget,  &
         530        1536 :                      history_budget_histfile_num_out = history_budget_histfile_num)
         531             :
         532             :    ! "irad_always" is number of time steps to execute radiation continuously from
         533             :    ! start of initial OR restart run
         534        1536 :    nstep = get_nstep()
         535        1536 :    if (irad_always > 0) then
         536           0 :       irad_always = irad_always + nstep
         537             :    end if
         538             :
         539        1536 :    if (docosp) call cospsimulator_intr_init()
         540             :
         541        4608 :    allocate(cosp_cnt(begchunk:endchunk), stat=istat)
         542        1536 :    call handle_allocate_error(istat, sub, 'cosp_cnt')
         543        1536 :    if (is_first_restart_step()) then
         544        3864 :       cosp_cnt(begchunk:endchunk) = cosp_cnt_init
         545             :    else
         546        3864 :       cosp_cnt(begchunk:endchunk) = 0
         547             :    end if
         548             :
         549             :    ! Add fields to history buffer
         550             :
         551             :    call addfld('TOT_CLD_VISTAU',  (/ 'lev' /), 'A',   '1',             &
         552             :                'Total gbx cloud extinction visible sw optical depth',  &
         553        3072 :                sampling_seq='rad_lwsw', flag_xyfill=.true.)
         554             :    call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A',  '1',              &
         555             :                'Total in-cloud extinction visible sw optical depth',   &
         556        3072 :                sampling_seq='rad_lwsw', flag_xyfill=.true.)
         557             :    call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A',  '1', &
         558             :                'Liquid in-cloud extinction visible sw optical depth',  &
         559        3072 :                sampling_seq='rad_lwsw', flag_xyfill=.true.)
         560             :    call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A',  '1',              &
         561             :                'Ice in-cloud extinction visible sw optical depth',     &
         562        3072 :                sampling_seq='rad_lwsw', flag_xyfill=.true.)
         563             :
         564        1536 :    if (cldfsnow_idx > 0) then
         565             :       call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1',           &
         566             :                   'Snow in-cloud extinction visible sw optical depth', &
         567        3072 :                   sampling_seq='rad_lwsw', flag_xyfill=.true.)
         568             :    end if
         569        1536 :    if (cldfgrau_idx > 0 .and. graupel_in_rad) then
         570             :       call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1',           &
         571             :                   'Graupel in-cloud extinction visible sw optical depth', &
         572           0 :                   sampling_seq='rad_lwsw', flag_xyfill=.true.)
         573             :    endif
         574             :
         575             :    ! get list of active radiation calls
         576        1536 :    call rad_cnst_get_call_list(active_calls)
         577             :
         578             :    ! Add shortwave radiation fields to history master field list.
         579             :
         580       18432 :    do icall = 0, N_DIAG
         581             :
         582       18432 :       if (active_calls(icall)) then
         583             :
         584             :          call addfld('SOLIN'//diag(icall),    horiz_only,   'A', 'W/m2', &
         585        1536 :                      'Solar insolation', sampling_seq='rad_lwsw')
         586             :          call addfld('QRS'//diag(icall),      (/ 'lev' /),  'A', 'K/s',  &
         587        3072 :                      'Solar heating rate', sampling_seq='rad_lwsw')
         588             :          call addfld('QRSC'//diag(icall),     (/ 'lev' /),  'A', 'K/s',  &
         589        3072 :                      'Clearsky solar heating rate', sampling_seq='rad_lwsw')
         590             :          call addfld('FSNT'//diag(icall),     horiz_only,   'A', 'W/m2', &
         591        1536 :                      'Net solar flux at top of model', sampling_seq='rad_lwsw')
         592             :          call addfld('FSNTC'//diag(icall),    horiz_only,   'A', 'W/m2', &
         593        1536 :                      'Clearsky net solar flux at top of model', sampling_seq='rad_lwsw')
         594             :          call addfld('FSNTOA'//diag(icall),   horiz_only,   'A', 'W/m2', &
         595        1536 :                      'Net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
         596             :          call addfld('FSNTOAC'//diag(icall),  horiz_only,   'A', 'W/m2', &
         597        1536 :                      'Clearsky net solar flux at top of atmosphere', sampling_seq='rad_lwsw')
         598             :          call addfld('SWCF'//diag(icall),     horiz_only,   'A', 'W/m2', &
         599        1536 :                      'Shortwave cloud forcing', sampling_seq='rad_lwsw')
         600             :          call addfld('FSUTOA'//diag(icall),   horiz_only,   'A', 'W/m2', &
         601        1536 :                      'Upwelling solar flux at top of atmosphere', sampling_seq='rad_lwsw')
         602             :          call addfld('FSNIRTOA'//diag(icall), horiz_only,   'A', 'W/m2', &
         603        1536 :                      'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
         604             :          call addfld('FSNRTOAC'//diag(icall), horiz_only,   'A', 'W/m2', &
         605        1536 :                       'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
         606             :          call addfld('FSNRTOAS'//diag(icall), horiz_only,   'A', 'W/m2', &
         607        1536 :                      'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw')
         608             :          call addfld('FSN200'//diag(icall),   horiz_only,   'A', 'W/m2', &
         609        1536 :                      'Net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
         610             :          call addfld('FSN200C'//diag(icall),  horiz_only,   'A', 'W/m2', &
         611        1536 :                      'Clearsky net shortwave flux at 200 mb', sampling_seq='rad_lwsw')
         612             :          call addfld('FSNR'//diag(icall),     horiz_only,   'A', 'W/m2', &
         613        1536 :                      'Net solar flux at tropopause', sampling_seq='rad_lwsw')
         614             :          call addfld('SOLL'//diag(icall),     horiz_only,   'A', 'W/m2', &
         615        1536 :                      'Solar downward near infrared direct  to surface', sampling_seq='rad_lwsw')
         616             :          call addfld('SOLS'//diag(icall),     horiz_only,   'A', 'W/m2', &
         617        1536 :                      'Solar downward visible direct  to surface', sampling_seq='rad_lwsw')
         618             :          call addfld('SOLLD'//diag(icall),    horiz_only,   'A', 'W/m2', &
         619        1536 :                      'Solar downward near infrared diffuse to surface', sampling_seq='rad_lwsw')
         620             :          call addfld('SOLSD'//diag(icall),    horiz_only,   'A', 'W/m2', &
         621        1536 :                      'Solar downward visible diffuse to surface', sampling_seq='rad_lwsw')
         622             :          call addfld('FSNS'//diag(icall),     horiz_only,   'A', 'W/m2', &
         623        1536 :                      'Net solar flux at surface', sampling_seq='rad_lwsw')
         624             :          call addfld('FSNSC'//diag(icall),    horiz_only,   'A', 'W/m2', &
         625        1536 :                      'Clearsky net solar flux at surface', sampling_seq='rad_lwsw')
         626             :          call addfld('FSDS'//diag(icall),     horiz_only,   'A', 'W/m2', &
         627        1536 :                      'Downwelling solar flux at surface', sampling_seq='rad_lwsw')
         628             :          call addfld('FSDSC'//diag(icall),    horiz_only,   'A', 'W/m2', &
         629        1536 :                      'Clearsky downwelling solar flux at surface', sampling_seq='rad_lwsw')
         630             :
         631             :          ! Fluxes on CAM grid
         632             :          call addfld('FUS'//diag(icall),      (/ 'ilev' /), 'I', 'W/m2', &
         633        3072 :                      'Shortwave upward flux', sampling_seq='rad_lwsw')
         634             :          call addfld('FDS'//diag(icall),      (/ 'ilev' /), 'I', 'W/m2', &
         635        3072 :                      'Shortwave downward flux', sampling_seq='rad_lwsw')
         636             :          call addfld('FUSC'//diag(icall),     (/ 'ilev' /), 'I', 'W/m2', &
         637        3072 :                      'Shortwave clear-sky upward flux', sampling_seq='rad_lwsw')
         638             :          call addfld('FDSC'//diag(icall),     (/ 'ilev' /), 'I', 'W/m2', &
         639        3072 :                      'Shortwave clear-sky downward flux', sampling_seq='rad_lwsw')
         640             :
         641        1536 :          if (history_amwg) then
         642        1536 :             call add_default('SOLIN'//diag(icall),   1, ' ')
         643        1536 :             call add_default('QRS'//diag(icall),     1, ' ')
         644        1536 :             call add_default('FSNT'//diag(icall),    1, ' ')
         645        1536 :             call add_default('FSNTC'//diag(icall),   1, ' ')
         646        1536 :             call add_default('FSNTOA'//diag(icall),  1, ' ')
         647        1536 :             call add_default('FSNTOAC'//diag(icall), 1, ' ')
         648        1536 :             call add_default('SWCF'//diag(icall),    1, ' ')
         649        1536 :             call add_default('FSNS'//diag(icall),    1, ' ')
         650        1536 :             call add_default('FSNSC'//diag(icall),   1, ' ')
         651        1536 :             call add_default('FSUTOA'//diag(icall),  1, ' ')
         652        1536 :             call add_default('FSDSC'//diag(icall),   1, ' ')
         653        1536 :             call add_default('FSDS'//diag(icall),    1, ' ')
         654             :          endif
         655             :
         656             :       end if
         657             :    end do
         658             :
         659        1536 :    if (scm_crm_mode) then
         660           0 :       call add_default('FUS     ', 1, ' ')
         661           0 :       call add_default('FUSC    ', 1, ' ')
         662           0 :       call add_default('FDS     ', 1, ' ')
         663           0 :       call add_default('FDSC    ', 1, ' ')
         664             :    endif
         665             :
         666             :    ! Add longwave radiation fields to history master field list.
         667             :
         668       18432 :    do icall = 0, N_DIAG
         669             :
         670       18432 :       if (active_calls(icall)) then
         671             :          call addfld('QRL'//diag(icall),     (/ 'lev' /), 'A', 'K/s',  &
         672        3072 :                      'Longwave heating rate', sampling_seq='rad_lwsw')
         673             :          call addfld('QRLC'//diag(icall),    (/ 'lev' /), 'A', 'K/s',  &
         674        3072 :                      'Clearsky longwave heating rate', sampling_seq='rad_lwsw')
         675             :          call addfld('FLNT'//diag(icall),    horiz_only,  'A', 'W/m2', &
         676        1536 :                      'Net longwave flux at top of model', sampling_seq='rad_lwsw')
         677             :          call addfld('FLNTC'//diag(icall),   horiz_only,  'A', 'W/m2', &
         678        1536 :                      'Clearsky net longwave flux at top of model', sampling_seq='rad_lwsw')
         679             :          call addfld('FLUT'//diag(icall),    horiz_only,  'A', 'W/m2', &
         680        1536 :                      'Upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
         681             :          call addfld('FLUTC'//diag(icall),   horiz_only,  'A', 'W/m2', &
         682        1536 :                      'Clearsky upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
         683             :          call addfld('LWCF'//diag(icall),    horiz_only,  'A', 'W/m2', &
         684        1536 :                      'Longwave cloud forcing', sampling_seq='rad_lwsw')
         685             :          call addfld('FLN200'//diag(icall),  horiz_only,  'A', 'W/m2', &
         686        1536 :                      'Net longwave flux at 200 mb', sampling_seq='rad_lwsw')
         687             :          call addfld('FLN200C'//diag(icall), horiz_only,  'A', 'W/m2', &
         688        1536 :                      'Clearsky net longwave flux at 200 mb', sampling_seq='rad_lwsw')
         689             :          call addfld('FLNR'//diag(icall),    horiz_only,  'A', 'W/m2', &
         690        1536 :                      'Net longwave flux at tropopause', sampling_seq='rad_lwsw')
         691             :          call addfld('FLNS'//diag(icall),    horiz_only,  'A', 'W/m2', &
         692        1536 :                      'Net longwave flux at surface', sampling_seq='rad_lwsw')
         693             :          call addfld('FLNSC'//diag(icall),   horiz_only,  'A', 'W/m2', &
         694        1536 :                      'Clearsky net longwave flux at surface', sampling_seq='rad_lwsw')
         695             :          call addfld('FLDS'//diag(icall),    horiz_only,  'A', 'W/m2', &
         696        1536 :                      'Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
         697             :          call addfld('FLDSC'//diag(icall),   horiz_only,  'A', 'W/m2', &
         698        1536 :                      'Clearsky Downwelling longwave flux at surface', sampling_seq='rad_lwsw')
         699             :
         700             :          ! Fluxes on CAM grid
         701             :          call addfld('FUL'//diag(icall),    (/ 'ilev' /), 'I', 'W/m2', &
         702        3072 :                      'Longwave upward flux', sampling_seq='rad_lwsw')
         703             :          call addfld('FDL'//diag(icall),    (/ 'ilev' /), 'I', 'W/m2', &
         704        3072 :                      'Longwave downward flux', sampling_seq='rad_lwsw')
         705             :          call addfld('FULC'//diag(icall),   (/ 'ilev' /), 'I', 'W/m2', &
         706        3072 :                      'Longwave clear-sky upward flux', sampling_seq='rad_lwsw')
         707             :          call addfld('FDLC'//diag(icall),   (/ 'ilev' /), 'I', 'W/m2', &
         708        3072 :                      'Longwave clear-sky downward flux', sampling_seq='rad_lwsw')
         709             :
         710        1536 :          if (history_amwg) then
         711        1536 :             call add_default('QRL'//diag(icall),   1, ' ')
         712        1536 :             call add_default('FLNT'//diag(icall),  1, ' ')
         713        1536 :             call add_default('FLNTC'//diag(icall), 1, ' ')
         714        1536 :             call add_default('FLUT'//diag(icall),  1, ' ')
         715        1536 :             call add_default('FLUTC'//diag(icall), 1, ' ')
         716        1536 :             call add_default('LWCF'//diag(icall),  1, ' ')
         717        1536 :             call add_default('FLNS'//diag(icall),  1, ' ')
         718        1536 :             call add_default('FLNSC'//diag(icall), 1, ' ')
         719        1536 :             call add_default('FLDS'//diag(icall),  1, ' ')
         720             :          end if
         721             :
         722             :       end if
         723             :    end do
         724             :
         725        3072 :    call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity')
         726             :
         727        1536 :    if (scm_crm_mode) then
         728           0 :       call add_default ('FUL     ', 1, ' ')
         729           0 :       call add_default ('FULC    ', 1, ' ')
         730           0 :       call add_default ('FDL     ', 1, ' ')
         731           0 :       call add_default ('FDLC    ', 1, ' ')
         732             :    endif
         733             :
         734             :    ! Heating rate needed for d(theta)/dt computation
         735        3072 :    call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation')
         736             :
         737        1536 :    if ( history_budget .and. history_budget_histfile_num > 1 ) then
         738           0 :       call add_default ('QRL     ', history_budget_histfile_num, ' ')
         739           0 :       call add_default ('QRS     ', history_budget_histfile_num, ' ')
         740             :    end if
         741             :
         742        1536 :    if (history_vdiag) then
         743           0 :       call add_default('FLUT', 2, ' ')
         744           0 :       call add_default('FLUT', 3, ' ')
         745             :    end if
         746             :
         747        3072 : end subroutine radiation_init
         748             :
         749             : !===============================================================================
         750             :
         751        1536 : subroutine radiation_define_restart(file)
         752             :
         753             :    ! define variables to be written to restart file
         754             :
         755             :    ! arguments
         756             :    type(file_desc_t), intent(inout) :: file
         757             :
         758             :    ! local variables
         759             :    integer :: ierr
         760             :    !----------------------------------------------------------------------------
         761             :
         762        1536 :    call pio_seterrorhandling(file, PIO_BCAST_ERROR)
         763             :
         764        1536 :    ierr = pio_def_var(file, 'nextsw_cday', pio_double, nextsw_cday_desc)
         765        1536 :    ierr = pio_put_att(file, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models')
         766        1536 :    if (docosp) then
         767           0 :       ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc)
         768             :    end if
         769             :
         770        1536 : end subroutine radiation_define_restart
         771             :
         772             : !===============================================================================
         773             :
         774        1536 : subroutine radiation_write_restart(file)
         775             :
         776             :    ! write variables to restart file
         777             :
         778             :    ! arguments
         779             :    type(file_desc_t), intent(inout) :: file
         780             :
         781             :    ! local variables
         782             :    integer :: ierr
         783             :    !----------------------------------------------------------------------------
         784        3072 :    ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /))
         785        1536 :    if (docosp) then
         786           0 :       ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/))
         787             :    end if
         788             :
         789        1536 : end subroutine radiation_write_restart
         790             :
         791             : !===============================================================================
         792             :
         793         768 : subroutine radiation_read_restart(file)
         794             :
         795             :    ! read variables from restart file
         796             :
         797             :    ! arguments
         798             :    type(file_desc_t), intent(inout) :: file
         799             :
         800             :    ! local variables
         801             :    integer :: ierr
         802             :    type(var_desc_t) :: vardesc
         803             :    integer :: err_handling
         804             :
         805             :    !----------------------------------------------------------------------------
         806         768 :    if (docosp) then
         807           0 :       call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
         808           0 :       ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc)
         809           0 :       call pio_seterrorhandling(File, err_handling)
         810           0 :       if (ierr /= PIO_NOERR) then
         811           0 :          cosp_cnt_init = 0
         812             :       else
         813           0 :          ierr = pio_get_var(File, vardesc, cosp_cnt_init)
         814             :       end if
         815             :    end if
         816             :
         817         768 :    ierr = pio_inq_varid(file, 'nextsw_cday', vardesc)
         818         768 :    ierr = pio_get_var(file, vardesc, nextsw_cday)
         819             :
         820             :
         821         768 : end subroutine radiation_read_restart
         822             :
         823             : !===============================================================================
         824             :
         825    62675424 : subroutine radiation_tend( &
         826             :    state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)
         827             :
         828             :    !-----------------------------------------------------------------------
         829             :    !
         830             :    ! CAM driver for radiation computation.
         831             :    !
         832             :    !-----------------------------------------------------------------------
         833             :
         834             :    ! Location/Orbital Parameters for cosine zenith angle
         835             :    use phys_grid,          only: get_rlat_all_p, get_rlon_all_p
         836             :    use cam_control_mod,    only: eccen, mvelpp, lambm0, obliqr
         837             :    use shr_orb_mod,        only: shr_orb_decl, shr_orb_cosz
         838             :
         839             :    use rrtmgp_inputs,      only: rrtmgp_set_state, rrtmgp_set_gases_lw, rrtmgp_set_cloud_lw, &
         840             :                                  rrtmgp_set_aer_lw, rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, &
         841             :                                  rrtmgp_set_aer_sw
         842             :
         843             :    ! RRTMGP drivers for flux calculations.
         844             :    use mo_rte_lw,          only: rte_lw
         845             :    use mo_rte_sw,          only: rte_sw
         846             :
         847             :    use radheat,            only: radheat_tend
         848             :
         849             :    use radiation_data,     only: rad_data_write
         850             :
         851             :    use interpolate_data,   only: vertinterp
         852             :    use tropopause,         only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
         853             :    use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps
         854             :
         855             :
         856             :    ! Arguments
         857             :    type(physics_state), intent(in), target :: state
         858             :    type(physics_ptend), intent(out)        :: ptend
         859             :    type(physics_buffer_desc), pointer      :: pbuf(:)
         860             :    type(cam_out_t),     intent(inout)      :: cam_out
         861             :    type(cam_in_t),      intent(in)         :: cam_in
         862             :    real(r8),            intent(out)        :: net_flx(pcols)
         863             :
         864             :    type(rad_out_t), target, optional, intent(out) :: rd_out
         865             :
         866             :
         867             :    ! Local variables
         868             :    type(rad_out_t), pointer :: rd  ! allow rd_out to be optional by allocating a local object
         869             :                                    ! if the argument is not present
         870             :    logical  :: write_output
         871             :
         872             :    integer  :: i, k, istat
         873             :    integer  :: lchnk, ncol
         874             :    logical  :: dosw, dolw
         875             :    integer  :: icall           ! loop index for climate/diagnostic radiation calls
         876             :
         877             :    real(r8) :: calday          ! current calendar day
         878             :    real(r8) :: delta           ! Solar declination angle  in radians
         879             :    real(r8) :: eccf            ! Earth orbit eccentricity factor
         880             :    real(r8) :: clat(pcols)     ! current latitudes(radians)
         881             :    real(r8) :: clon(pcols)     ! current longitudes(radians)
         882             :    real(r8) :: coszrs(pcols)   ! Cosine solar zenith angle
         883             :
         884             :    ! Gathered indices of day and night columns
         885             :    !  chunk_column_index = IdxDay(daylight_column_index)
         886             :    integer :: Nday           ! Number of daylight columns
         887             :    integer :: Nnite          ! Number of night columns
         888             :    integer :: IdxDay(pcols)  ! chunk indices of daylight columns
         889             :    integer :: IdxNite(pcols) ! chunk indices of night columns
         890             :
         891             :    integer :: itim_old
         892             :
         893     1492272 :    real(r8), pointer :: cld(:,:)      ! cloud fraction
         894     1492272 :    real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
         895     1492272 :    real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
         896             :    real(r8)          :: cldfprime(pcols,pver)   ! combined cloud fraction
         897     1492272 :    real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate
         898     1492272 :    real(r8), pointer :: qrl(:,:) ! longwave  radiative heating rate
         899     1492272 :    real(r8), pointer :: fsds(:)  ! Surface solar down flux
         900     1492272 :    real(r8), pointer :: fsns(:)  ! Surface solar absorbed flux
         901     1492272 :    real(r8), pointer :: fsnt(:)  ! Net column abs solar flux at model top
         902     1492272 :    real(r8), pointer :: flns(:)  ! Srf longwave cooling (up-down) flux
         903     1492272 :    real(r8), pointer :: flnt(:)  ! Net outgoing lw flux at model top
         904             :
         905             :    real(r8), pointer, dimension(:,:,:) :: su => NULL()  ! shortwave spectral flux up
         906             :    real(r8), pointer, dimension(:,:,:) :: sd => NULL()  ! shortwave spectral flux down
         907             :    real(r8), pointer, dimension(:,:,:) :: lu => NULL()  ! longwave  spectral flux up
         908             :    real(r8), pointer, dimension(:,:,:) :: ld => NULL()  ! longwave  spectral flux down
         909             :
         910             :    ! tropopause diagnostic
         911             :    integer  :: troplev(pcols)
         912             :    real(r8) :: p_trop(pcols)
         913             :
         914             :    ! state data passed to radiation calc
         915     1492272 :    real(r8), allocatable :: t_sfc(:)
         916     1492272 :    real(r8), allocatable :: emis_sfc(:,:)
         917     1492272 :    real(r8), allocatable :: t_rad(:,:)
         918     1492272 :    real(r8), allocatable :: pmid_rad(:,:)
         919     1492272 :    real(r8), allocatable :: pint_rad(:,:)
         920     1492272 :    real(r8), allocatable :: t_day(:,:)
         921     1492272 :    real(r8), allocatable :: pmid_day(:,:)
         922     1492272 :    real(r8), allocatable :: pint_day(:,:)
         923     1492272 :    real(r8), allocatable :: coszrs_day(:)
         924     1492272 :    real(r8), allocatable :: alb_dir(:,:)
         925     1492272 :    real(r8), allocatable :: alb_dif(:,:)
         926             :
         927             :    ! in-cloud optical depths for COSP
         928             :    real(r8) :: cld_tau_cloudsim(pcols,pver)    ! liq + ice
         929             :    real(r8) :: snow_tau_cloudsim(pcols,pver)   ! snow
         930             :    real(r8) :: grau_tau_cloudsim(pcols,pver)   ! graupel
         931             :    real(r8) :: cld_lw_abs_cloudsim(pcols,pver) ! liq + ice
         932             :    real(r8) :: snow_lw_abs_cloudsim(pcols,pver)! snow
         933             :    real(r8) :: grau_lw_abs_cloudsim(pcols,pver)! graupel
         934             :
         935             :    ! Set vertical indexing in RRTMGP to be the same as CAM (top to bottom).
         936             :    logical, parameter :: top_at_1 = .true.
         937             :
         938             :    ! TOA solar flux on RRTMGP g-points
         939     1492272 :    real(r8), allocatable :: toa_flux(:,:)
         940             :    ! TSI from RRTMGP data (from sum over g-point representation)
         941             :    real(r8) :: tsi_ref
         942             :
         943             :    ! Planck sources for LW.
         944     1492272 :    type(ty_source_func_lw) :: sources_lw
         945             :
         946             :    ! Gas volume mixing ratios.  Use separate objects for LW and SW because SW only does
         947             :    ! calculations for daylight columns.
         948             :    ! These objects have a final method which deallocates the internal memory when they
         949             :    ! go out of scope (i.e., when radiation_tend returns), so no need for explicit deallocation.
         950     1492272 :    type(ty_gas_concs) :: gas_concs_lw
         951     1492272 :    type(ty_gas_concs) :: gas_concs_sw
         952             :
         953             :    ! Atmosphere optics.  This object is initialized with gas optics, then is incremented
         954             :    ! by the aerosol optics for the clear-sky radiative flux calculations, and then
         955             :    ! incremented again by the cloud optics for the all-sky radiative flux calculations.
         956     1492272 :    type(ty_optical_props_1scl) :: atm_optics_lw
         957     1492272 :    type(ty_optical_props_2str) :: atm_optics_sw
         958             :
         959             :    ! Cloud optical properties objects (McICA sampling of cloud optical properties).
         960     1492272 :    type(ty_optical_props_1scl) :: cloud_lw
         961     1492272 :    type(ty_optical_props_2str) :: cloud_sw
         962             :
         963             :    ! Aerosol optical properties objects.
         964     1492272 :    type(ty_optical_props_1scl) :: aer_lw
         965     1492272 :    type(ty_optical_props_2str) :: aer_sw
         966             :
         967             :    ! Flux objects contain all fluxes computed by RRTMGP.
         968             :    ! SW allsky fluxes always include spectrally resolved fluxes needed for surface models.
         969             :    type(ty_fluxes_byband) :: fsw
         970             :    ! LW allsky fluxes only need spectrally resolved fluxes when spectralflux=.true.
         971             :    type(ty_fluxes_byband) :: flw
         972             :    ! Only broadband fluxes needed for clear sky (diagnostics).
         973             :    type(ty_fluxes_broadband) :: fswc, flwc
         974             :
         975             :    ! Arrays for output diagnostics on CAM grid.
         976             :    real(r8) :: fns(pcols,pverp)     ! net shortwave flux
         977             :    real(r8) :: fcns(pcols,pverp)    ! net clear-sky shortwave flux
         978             :    real(r8) :: fnl(pcols,pverp)     ! net longwave flux
         979             :    real(r8) :: fcnl(pcols,pverp)    ! net clear-sky longwave flux
         980             :
         981             :    ! for COSP
         982             :    real(r8) :: emis(pcols,pver)        ! Cloud longwave emissivity
         983             :    real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau
         984             :    real(r8) :: gb_snow_lw(pcols,pver)  ! grid-box mean LW snow optical depth
         985             :
         986             :    real(r8) :: ftem(pcols,pver)        ! Temporary workspace for outfld variables
         987             :
         988             :    character(len=128) :: errmsg
         989             :    character(len=*), parameter :: sub = 'radiation_tend'
         990             :    !--------------------------------------------------------------------------------------
         991             :
         992     1492272 :    lchnk = state%lchnk
         993     1492272 :    ncol = state%ncol
         994             :
         995     1492272 :    if (present(rd_out)) then
         996           0 :       rd => rd_out
         997           0 :       write_output = .false.
         998             :    else
         999     1492272 :       allocate(rd, stat=istat)
        1000     1492272 :       call handle_allocate_error(istat, sub, 'rd')
        1001     1492272 :       write_output = .true.
        1002             :    end if
        1003             :
        1004     1492272 :    dosw = radiation_do('sw', get_nstep())      ! do shortwave radiation calc this timestep?
        1005     1492272 :    dolw = radiation_do('lw', get_nstep())      ! do longwave radiation calc this timestep?
        1006             :
        1007             :    ! Cosine solar zenith angle for current time step
        1008     1492272 :    calday = get_curr_calday()
        1009     1492272 :    call get_rlat_all_p(lchnk, ncol, clat)
        1010     1492272 :    call get_rlon_all_p(lchnk, ncol, clon)
        1011             :
        1012             :    call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, &
        1013     1492272 :                      delta, eccf)
        1014             :
        1015     1492272 :    if (use_rad_uniform_angle) then
        1016           0 :       do i = 1, ncol
        1017           0 :          coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, &
        1018           0 :                                   uniform_angle=rad_uniform_angle)
        1019             :       end do
        1020             :    else
        1021    24917472 :       do i = 1, ncol
        1022             :          ! if dt_avg /= 0, it triggers using avg coszrs
        1023    24917472 :          coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg)
        1024             :       end do
        1025             :    end if
        1026             :
        1027             :    ! Gather night/day column indices.
        1028     1492272 :    Nday = 0
        1029     1492272 :    Nnite = 0
        1030    24917472 :    do i = 1, ncol
        1031    24917472 :       if ( coszrs(i) > 0.0_r8 ) then
        1032    11712600 :          Nday = Nday + 1
        1033    11712600 :          IdxDay(Nday) = i
        1034             :       else
        1035    11712600 :          Nnite = Nnite + 1
        1036    11712600 :          IdxNite(Nnite) = i
        1037             :       end if
        1038             :    end do
        1039             :
        1040             :    ! Associate pointers to physics buffer fields
        1041     1492272 :    itim_old = pbuf_old_tim_idx()
        1042     1492272 :    nullify(cldfsnow)
        1043     1492272 :    if (cldfsnow_idx > 0) then
        1044     5969088 :       call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
        1045             :    end if
        1046     1492272 :    nullify(cldfgrau)
        1047     1492272 :    if (cldfgrau_idx > 0 .and. graupel_in_rad) then
        1048           0 :       call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
        1049             :    endif
        1050             :
        1051     5969088 :    call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
        1052             :
        1053     1492272 :    call pbuf_get_field(pbuf, qrs_idx, qrs)
        1054     1492272 :    call pbuf_get_field(pbuf, qrl_idx, qrl)
        1055             :
        1056     1492272 :    call pbuf_get_field(pbuf, fsns_idx, fsns)
        1057     1492272 :    call pbuf_get_field(pbuf, fsnt_idx, fsnt)
        1058     1492272 :    call pbuf_get_field(pbuf, fsds_idx, fsds)
        1059     1492272 :    call pbuf_get_field(pbuf, flns_idx, flns)
        1060     1492272 :    call pbuf_get_field(pbuf, flnt_idx, flnt)
        1061             :
        1062     1492272 :    if (spectralflux) then
        1063           0 :       call pbuf_get_field(pbuf, su_idx, su)
        1064           0 :       call pbuf_get_field(pbuf, sd_idx, sd)
        1065           0 :       call pbuf_get_field(pbuf, lu_idx, lu)
        1066           0 :       call pbuf_get_field(pbuf, ld_idx, ld)
        1067             :    end if
        1068             :
        1069             :    ! Allocate the flux arrays and init to zero.
        1070     1492272 :    call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fsw, do_direct=.true.)
        1071     1492272 :    call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fswc, do_direct=.true.)
        1072     1492272 :    call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flw)
        1073     1492272 :    call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flwc)
        1074             :
        1075             :    !  For CRM, make cloud equal to input observations:
        1076     1492272 :    if (scm_crm_mode .and. have_cld) then
        1077           0 :       do k = 1, pver
        1078           0 :          cld(:ncol,k)= cldobs(k)
        1079             :       end do
        1080             :    end if
        1081             :
        1082             :    ! Find tropopause height if needed for diagnostic output
        1083     1492272 :    if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
        1084             :       !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
        1085           0 :       troplev(:) = 0
        1086           0 :       p_trop(:) = 0._r8
        1087             :       !REMOVECAM_END
        1088             :       call tropopause_find_cam(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, &
        1089           0 :                            backup=TROP_ALG_CLIMATE)
        1090             :    end if
        1091             :
        1092             :    ! Get time of next radiation calculation - albedos will need to be
        1093             :    ! calculated by each surface model at this time
        1094     1492272 :    nextsw_cday = radiation_nextsw_cday()
        1095             :
        1096     1492272 :    if (dosw .or. dolw) then
        1097             :
        1098             :       allocate( &
        1099             :          t_sfc(ncol), emis_sfc(nlwbands,ncol), toa_flux(nday,nswgpts),     &
        1100             :          t_rad(ncol,nlay), pmid_rad(ncol,nlay), pint_rad(ncol,nlay+1),     &
        1101             :          t_day(nday,nlay), pmid_day(nday,nlay), pint_day(nday,nlay+1),     &
        1102             :          coszrs_day(nday), alb_dir(nswbands,nday), alb_dif(nswbands,nday), &
        1103    19139833 :          stat=istat)
        1104      746136 :       call handle_allocate_error(istat, sub, 't_sfc,..,alb_dif')
        1105             :
        1106             :       ! Prepares state variables, daylit columns, albedos for RRTMGP
        1107             :       call rrtmgp_set_state( &
        1108             :          state, cam_in, ncol, nlay, nday,             &
        1109             :          idxday, coszrs, kdist_sw, t_sfc, emis_sfc,   &
        1110             :          t_rad, pmid_rad, pint_rad, t_day, pmid_day,  &
        1111      746136 :          pint_day, coszrs_day, alb_dir, alb_dif)
        1112             :
        1113             :       ! Output the mass per layer, and total column burdens for gas and aerosol
        1114             :       ! constituents in the climate list.
        1115      746136 :       call rad_cnst_out(0, state, pbuf)
        1116             :
        1117             :       ! Modified cloud fraction accounts for radiatively active snow and/or graupel
        1118      746136 :       call modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
        1119             :
        1120             :       !========================!
        1121             :       ! SHORTWAVE calculations !
        1122             :       !========================!
        1123             :
        1124      746136 :       if (dosw) then
        1125             :
        1126             :          ! Set cloud optical properties in cloud_sw object.
        1127             :          call rrtmgp_set_cloud_sw( &
        1128             :             state, pbuf, nlay, nday, idxday,                              &
        1129             :             nnite, idxnite, pmid_day, cld, cldfsnow,                      &
        1130             :             cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw,      &
        1131             :             rd%tot_cld_vistau, rd%tot_icld_vistau, rd%liq_icld_vistau,    &
        1132             :             rd%ice_icld_vistau, rd%snow_icld_vistau, rd%grau_icld_vistau, &
        1133      746136 :             cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim )
        1134             :
        1135      746136 :          if (write_output) then
        1136      746136 :             call radiation_output_cld(lchnk, rd)
        1137             :          end if
        1138             :
        1139             :          ! If no daylight columns, can't create empty RRTMGP objects
        1140      746136 :          if (nday > 0) then
        1141             :
        1142             :             ! Initialize object for gas concentrations.
        1143      389263 :             errmsg = gas_concs_sw%init(gaslist_lc)
        1144      389263 :             call stop_on_err(errmsg, sub, 'gas_concs_sw%init')
        1145             :
        1146             :             ! Initialize object for combined gas + aerosol + cloud optics.
        1147             :             ! Allocates arrays for properties represented on g-points.
        1148      389263 :             errmsg = atm_optics_sw%alloc_2str(nday, nlay, kdist_sw)
        1149      389263 :             call stop_on_err(errmsg, sub, 'atm_optics_sw%alloc_2str')
        1150             :
        1151             :             ! Initialize object for SW aerosol optics.  Allocates arrays
        1152             :             ! for properties represented by band.
        1153      389263 :             errmsg = aer_sw%alloc_2str(nday, nlay, kdist_sw%get_band_lims_wavenumber())
        1154      389263 :             call stop_on_err(errmsg, sub, 'aer_sw%alloc_2str')
        1155             :
        1156             :          end if
        1157             :
        1158             :          ! The climate (icall==0) calculation must occur last.
        1159     8953632 :          do icall = N_DIAG, 0, -1
        1160     8953632 :             if (active_calls(icall)) then
        1161             :
        1162      746136 :                if (nday > 0) then
        1163             :
        1164             :                   ! Set gas volume mixing ratios for this call in gas_concs_sw.
        1165             :                   call rrtmgp_set_gases_sw( &
        1166             :                      icall, state, pbuf, nlay, nday, &
        1167      389263 :                      idxday, gas_concs_sw)
        1168             :
        1169             :                   ! Compute the gas optics (stored in atm_optics_sw).
        1170             :                   ! toa_flux is the reference solar source from RRTMGP data.
        1171             :                   errmsg = kdist_sw%gas_optics( &
        1172             :                      pmid_day, pint_day, t_day, gas_concs_sw, atm_optics_sw, &
        1173      389263 :                      toa_flux)
        1174      389263 :                   call stop_on_err(errmsg, sub, 'kdist_sw%gas_optics')
        1175             :
        1176             :                   ! Scale the solar source
        1177    43986719 :                   tsi_ref = sum(toa_flux(1,:))
        1178   699892319 :                   toa_flux = toa_flux * sol_tsi * eccf / tsi_ref
        1179             :
        1180             :                end if
        1181             :
        1182             :                ! Set SW aerosol optical properties in the aer_sw object.
        1183             :                ! This call made even when no daylight columns because it does some
        1184             :                ! diagnostic aerosol output.
        1185             :                call rrtmgp_set_aer_sw( &
        1186      746136 :                   icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
        1187             :
        1188      746136 :                if (nday > 0) then
        1189             :
        1190             :                   ! Increment the gas optics (in atm_optics_sw) by the aerosol optics in aer_sw.
        1191      389263 :                   errmsg = aer_sw%increment(atm_optics_sw)
        1192      389263 :                   call stop_on_err(errmsg, sub, 'aer_sw%increment')
        1193             :
        1194             :                   ! Compute clear-sky fluxes.
        1195             :                   errmsg = rte_sw(&
        1196             :                      atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
        1197      389263 :                      alb_dir, alb_dif, fswc)
        1198      389263 :                   call stop_on_err(errmsg, sub, 'clear-sky rte_sw')
        1199             :
        1200             :                   ! Increment the aerosol+gas optics (in atm_optics_sw) by the cloud optics in cloud_sw.
        1201      389263 :                   errmsg = cloud_sw%increment(atm_optics_sw)
        1202      389263 :                   call stop_on_err(errmsg, sub, 'cloud_sw%increment')
        1203             :
        1204             :                   ! Compute all-sky fluxes.
        1205             :                   errmsg = rte_sw(&
        1206             :                      atm_optics_sw, top_at_1, coszrs_day, toa_flux, &
        1207      389263 :                      alb_dir, alb_dif, fsw)
        1208      389263 :                   call stop_on_err(errmsg, sub, 'all-sky rte_sw')
        1209             :
        1210             :                end if
        1211             :
        1212             :                ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
        1213      746136 :                call set_sw_diags()
        1214             :
        1215      746136 :                if (write_output) then
        1216      746136 :                   call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
        1217             :                end if
        1218             :
        1219             :             end if ! (active_calls(icall))
        1220             :          end do    ! loop over diagnostic calcs (icall)
        1221             :
        1222             :       else
        1223             :          ! SW calc not done.  pbuf carries Q*dp across timesteps.
        1224             :          ! Convert to Q before calling radheat_tend.
        1225           0 :          qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
        1226             :       end if  ! if (dosw)
        1227             :
        1228             :       !=======================!
        1229             :       ! LONGWAVE calculations !
        1230             :       !=======================!
        1231             :
        1232      746136 :       if (dolw) then
        1233             :
        1234             :          ! Initialize object for Planck sources.
        1235      746136 :          errmsg = sources_lw%alloc(ncol, nlay, kdist_lw)
        1236      746136 :          call stop_on_err(errmsg, sub, 'sources_lw%alloc')
        1237             :
        1238             :          ! Set cloud optical properties in cloud_lw object.
        1239             :          call rrtmgp_set_cloud_lw( &
        1240             :             state, pbuf, ncol, nlay, nlaycam, &
        1241             :             cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
        1242      746136 :             kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
        1243             :
        1244             :          ! Initialize object for gas concentrations
        1245      746136 :          errmsg = gas_concs_lw%init(gaslist_lc)
        1246      746136 :          call stop_on_err(errmsg, sub, 'gas_concs_lw%init')
        1247             :
        1248             :          ! Initialize object for combined gas + aerosol + cloud optics.
        1249      746136 :          errmsg = atm_optics_lw%alloc_1scl(ncol, nlay, kdist_lw)
        1250      746136 :          call stop_on_err(errmsg, sub, 'atm_optics_lw%alloc_1scl')
        1251             :
        1252             :          ! Initialize object for LW aerosol optics.
        1253      746136 :          errmsg = aer_lw%alloc_1scl(ncol, nlay, kdist_lw%get_band_lims_wavenumber())
        1254      746136 :          call stop_on_err(errmsg, sub, 'aer_lw%alloc_1scl')
        1255             :
        1256             :          ! The climate (icall==0) calculation must occur last.
        1257     8953632 :          do icall = N_DIAG, 0, -1
        1258             :
        1259     8953632 :             if (active_calls(icall)) then
        1260             :
        1261             :                ! Set gas volume mixing ratios for this call in gas_concs_lw.
        1262      746136 :                call rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs_lw)
        1263             :
        1264             :                ! Compute the gas optics and Planck sources.
        1265             :                errmsg = kdist_lw%gas_optics( &
        1266             :                   pmid_rad, pint_rad, t_rad, t_sfc, gas_concs_lw, &
        1267      746136 :                   atm_optics_lw, sources_lw)
        1268      746136 :                call stop_on_err(errmsg, sub, 'kdist_lw%gas_optics')
        1269             :
        1270             :                ! Set LW aerosol optical properties in the aer_lw object.
        1271      746136 :                call rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
        1272             :
        1273             :                ! Increment the gas optics by the aerosol optics.
        1274      746136 :                errmsg = aer_lw%increment(atm_optics_lw)
        1275      746136 :                call stop_on_err(errmsg, sub, 'aer_lw%increment')
        1276             :
        1277             :                ! Compute clear-sky LW fluxes
        1278      746136 :                errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flwc)
        1279      746136 :                call stop_on_err(errmsg, sub, 'clear-sky rte_lw')
        1280             :
        1281             :                ! Increment the gas+aerosol optics by the cloud optics.
        1282      746136 :                errmsg = cloud_lw%increment(atm_optics_lw)
        1283      746136 :                call stop_on_err(errmsg, sub, 'cloud_lw%increment')
        1284             :
        1285             :                ! Compute all-sky LW fluxes
        1286      746136 :                errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flw)
        1287      746136 :                call stop_on_err(errmsg, sub, 'all-sky rte_lw')
        1288             :
        1289             :                ! Transform RRTMGP outputs to CAM outputs and compute heating rates.
        1290      746136 :                call set_lw_diags()
        1291             :
        1292      746136 :                if (write_output) then
        1293      746136 :                   call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
        1294             :                end if
        1295             :
        1296             :             end if ! (active_calls(icall))
        1297             :          end do    ! loop over diagnostic calcs (icall)
        1298             :
        1299             :       else
        1300             :          ! LW calc not done.  pbuf carries Q*dp across timesteps.
        1301             :          ! Convert to Q before calling radheat_tend.
        1302           0 :         qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
        1303             :       end if  ! if (dolw)
        1304             :
        1305           0 :       deallocate( &
        1306           0 :          t_sfc, emis_sfc, toa_flux, t_rad, pmid_rad, pint_rad,  &
        1307      746136 :          t_day, pmid_day, pint_day, coszrs_day, alb_dir, alb_dif)
        1308             :
        1309             :       !================!
        1310             :       ! COSP simulator !
        1311             :       !================!
        1312             :
        1313      746136 :       if (docosp) then
        1314             :
        1315           0 :          emis(:,:) = 0._r8
        1316           0 :          emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs_cloudsim(:ncol,:))
        1317           0 :          call outfld('EMIS', emis, pcols, lchnk)
        1318             :
        1319             :          ! compute grid-box mean SW and LW snow optical depth for use by COSP
        1320           0 :          gb_snow_tau(:,:) = 0._r8
        1321           0 :          gb_snow_lw(:,:)  = 0._r8
        1322           0 :          if (cldfsnow_idx > 0) then
        1323           0 :             do i = 1, ncol
        1324           0 :                do k = 1, pver
        1325           0 :                   if (cldfsnow(i,k) > 0._r8) then
        1326             :
        1327             :                      ! Add graupel to snow tau for cosp
        1328           0 :                      if (cldfgrau_idx > 0 .and. graupel_in_rad) then
        1329           0 :                         gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k) + &
        1330           0 :                                            grau_tau_cloudsim(i,k)*cldfgrau(i,k)
        1331             :                         gb_snow_lw(i,k)  = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k) + &
        1332           0 :                                            grau_lw_abs_cloudsim(i,k)*cldfgrau(i,k)
        1333             :                      else
        1334           0 :                         gb_snow_tau(i,k) = snow_tau_cloudsim(i,k)*cldfsnow(i,k)
        1335           0 :                         gb_snow_lw(i,k)  = snow_lw_abs_cloudsim(i,k)*cldfsnow(i,k)
        1336             :                      end if
        1337             :                   end if
        1338             :                end do
        1339             :             end do
        1340             :          end if
        1341             :
        1342             :          ! advance counter for this timestep (chunk dimension required for thread safety)
        1343           0 :          cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1
        1344             :
        1345             :          ! if counter is the same as cosp_nradsteps, run cosp and reset counter
        1346           0 :          if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then
        1347             :
        1348             :             ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave
        1349             :             !       optical depths are passed.
        1350             :             call cospsimulator_intr_run( &
        1351             :                state,  pbuf, cam_in, emis, coszrs,                     &
        1352             :                cld_swtau_in=cld_tau_cloudsim, snow_tau_in=gb_snow_tau, &
        1353           0 :                snow_emis_in=gb_snow_lw)
        1354           0 :             cosp_cnt(lchnk) = 0
        1355             :          end if
        1356             :       end if   ! docosp
        1357             :
        1358             :    else
        1359             :       ! Radiative flux calculations not done.  The quantity Q*dp is carried by the
        1360             :       ! physics buffer across timesteps.  It must be converted to Q (dry static energy
        1361             :       ! tendency) before being passed to radheat_tend.
        1362  1159408584 :       qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:)
        1363  1159408584 :       qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:)
        1364             :
        1365             :    end if   ! if (dosw .or. dolw) then
        1366             :
        1367             :    ! Output for PORT: Parallel Offline Radiative Transport
        1368     1492272 :    call rad_data_write(pbuf, state, cam_in, coszrs)
        1369             :
        1370             :    ! Compute net radiative heating tendency.  Note that the WACCM version
        1371             :    ! of radheat_tend merges upper atmosphere heating rates with those calculated
        1372             :    ! by RRTMGP.
        1373             :    call radheat_tend(state, pbuf,  ptend, qrl, qrs, fsns, &
        1374     1492272 :                      fsnt, flns, flnt, cam_in%asdir, net_flx)
        1375             :
        1376     1492272 :    if (write_output) then
        1377             :       ! Compute heating rate for dtheta/dt
        1378   140273568 :       do k = 1, pver
        1379  2318817168 :          do i = 1, ncol
        1380  2317324896 :             ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
        1381             :          end do
        1382             :       end do
        1383     1492272 :       call outfld('HR', ftem, pcols, lchnk)
        1384             :    end if
        1385             :
        1386             :    ! The radiative heating rates are carried in the physics buffer across timesteps
        1387             :    ! as Q*dp (for energy conservation).
        1388  2318817168 :    qrs(:ncol,:) = qrs(:ncol,:) * state%pdel(:ncol,:)
        1389  2318817168 :    qrl(:ncol,:) = qrl(:ncol,:) * state%pdel(:ncol,:)
        1390             :
        1391     1492272 :    if (.not. present(rd_out)) then
        1392     1492272 :       deallocate(rd)
        1393             :    end if
        1394     1492272 :    call free_optics_sw(atm_optics_sw)
        1395     1492272 :    call free_optics_sw(cloud_sw)
        1396     1492272 :    call free_optics_sw(aer_sw)
        1397     1492272 :    call free_fluxes(fsw)
        1398     1492272 :    call free_fluxes(fswc)
        1399             :
        1400     1492272 :    call sources_lw%finalize()
        1401     1492272 :    call free_optics_lw(cloud_lw)
        1402     1492272 :    call free_optics_lw(aer_lw)
        1403     1492272 :    call free_fluxes(flw)
        1404     2984544 :    call free_fluxes(flwc)
        1405             :
        1406             :    !-------------------------------------------------------------------------------
        1407             :    contains
        1408             :    !-------------------------------------------------------------------------------
        1409             :
        1410      746136 :    subroutine set_sw_diags()
        1411             :
        1412             :       ! Transform RRTMGP output for CAM and compute heating rates.
        1413             :       ! SW fluxes from RRTMGP are on daylight columns only, so expand to
        1414             :       ! full chunks for output to CAM history.
        1415             :
        1416             :       integer :: i
        1417             :       real(r8), dimension(size(fsw%bnd_flux_dn,1), &
        1418             :                           size(fsw%bnd_flux_dn,2), &
        1419     1492272 :                           size(fsw%bnd_flux_dn,3)) :: flux_dn_diffuse
        1420             :       !-------------------------------------------------------------------------
        1421             :
        1422             :       ! Initialize to provide 0.0 values for night columns.
        1423      746136 :       fns               = 0._r8 ! net sw flux
        1424      746136 :       fcns              = 0._r8 ! net sw clearsky flux
        1425    12684312 :       fsds              = 0._r8 ! downward sw flux at surface
        1426    12684312 :       rd%fsdsc          = 0._r8 ! downward sw clearsky flux at surface
        1427    12684312 :       rd%fsutoa         = 0._r8 ! upward sw flux at TOA
        1428    12684312 :       rd%fsntoa         = 0._r8 ! net sw at TOA
        1429    12684312 :       rd%fsntoac        = 0._r8 ! net sw clearsky flux at TOA
        1430    12684312 :       rd%solin          = 0._r8 ! solar irradiance at TOA
        1431  1193071464 :       rd%flux_sw_up     = 0._r8
        1432  1193071464 :       rd%flux_sw_dn     = 0._r8
        1433  1193071464 :       rd%flux_sw_clr_up = 0._r8
        1434  1193071464 :       rd%flux_sw_clr_dn = 0._r8
        1435             :
        1436  1180387152 :       qrs      = 0._r8
        1437    12684312 :       fsns     = 0._r8
        1438    12684312 :       fsnt     = 0._r8
        1439  1180387152 :       rd%qrsc  = 0._r8
        1440    12684312 :       rd%fsnsc = 0._r8
        1441    12684312 :       rd%fsntc = 0._r8
        1442             :
        1443     6602436 :       do i = 1, nday
        1444   556348500 :          fns(idxday(i),ktopcam:)  = fsw%flux_net(i, ktoprad:)
        1445   556348500 :          fcns(idxday(i),ktopcam:) = fswc%flux_net(i,ktoprad:)
        1446     5856300 :          fsds(idxday(i))          = fsw%flux_dn(i, nlay+1)
        1447     5856300 :          rd%fsdsc(idxday(i))      = fswc%flux_dn(i, nlay+1)
        1448     5856300 :          rd%fsutoa(idxday(i))     = fsw%flux_up(i, 1)
        1449     5856300 :          rd%fsntoa(idxday(i))     = fsw%flux_net(i, 1)
        1450     5856300 :          rd%fsntoac(idxday(i))    = fswc%flux_net(i, 1)
        1451     5856300 :          rd%solin(idxday(i))      = fswc%flux_dn(i, 1)
        1452   556348500 :          rd%flux_sw_up(idxday(i),ktopcam:)     = fsw%flux_up(i,ktoprad:)
        1453   556348500 :          rd%flux_sw_dn(idxday(i),ktopcam:)     = fsw%flux_dn(i,ktoprad:)
        1454   556348500 :          rd%flux_sw_clr_up(idxday(i),ktopcam:) = fswc%flux_up(i,ktoprad:)
        1455   557094636 :          rd%flux_sw_clr_dn(idxday(i),ktopcam:) = fswc%flux_dn(i,ktoprad:)
        1456             :       end do
        1457             :
        1458             :       ! Compute heating rate as a dry static energy tendency.
        1459      746136 :       call heating_rate('SW', ncol, fns, qrs)
        1460      746136 :       call heating_rate('SW', ncol, fcns, rd%qrsc)
        1461             :
        1462    12458736 :       fsns(:ncol)     = fns(:ncol,pverp)    ! net sw flux at surface
        1463    12458736 :       fsnt(:ncol)     = fns(:ncol,ktopcam)  ! net sw flux at top-of-model (w/o extra layer)
        1464    13204872 :       rd%fsnsc(:ncol) = fcns(:ncol,pverp)   ! net sw clearsky flux at surface
        1465    13204872 :       rd%fsntc(:ncol) = fcns(:ncol,ktopcam) ! net sw clearsky flux at top
        1466             :
        1467    12458736 :       cam_out%netsw(:ncol) = fsns(:ncol)
        1468             :
        1469             :       ! Output fluxes at 200 mb
        1470      746136 :       call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns,  rd%fsn200)
        1471      746136 :       call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
        1472      746136 :       if (hist_fld_active('FSNR')) then
        1473           0 :          do i = 1,ncol
        1474           0 :             call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
        1475             :          end do
        1476             :       end if
        1477             :
        1478      746136 :       if (spectralflux) then
        1479           0 :          su  = 0._r8
        1480           0 :          sd  = 0._r8
        1481           0 :          do i = 1, nday
        1482           0 :             su(idxday(i),ktopcam:,:) = fsw%bnd_flux_up(i,ktoprad:,:)
        1483           0 :             sd(idxday(i),ktopcam:,:) = fsw%bnd_flux_dn(i,ktoprad:,:)
        1484             :          end do
        1485             :       end if
        1486             :
        1487             :       ! Export surface fluxes
        1488             :       ! sols(pcols)      Direct solar rad on surface (< 0.7)
        1489             :       ! soll(pcols)      Direct solar rad on surface (>= 0.7)
        1490             :       ! RRTMGP: Near-IR bands (1-10), 820-16000 cm-1, 0.625-12.195 microns
        1491             :       ! Put half of band 10 in each of the UV/visible and near-IR values,
        1492             :       ! since this band straddles 0.7 microns:
        1493             :       ! UV/visible bands 10-13, 16000-50000 cm-1, 0.200-0.625 micron
        1494             :
        1495             :       ! reset fluxes
        1496    12684312 :       cam_out%sols  = 0.0_r8
        1497    12684312 :       cam_out%soll  = 0.0_r8
        1498    12684312 :       cam_out%solsd = 0.0_r8
        1499    12684312 :       cam_out%solld = 0.0_r8
        1500             :
        1501             :       ! Calculate diffuse flux from total and direct
        1502  8792431920 :       flux_dn_diffuse = fsw%bnd_flux_dn - fsw%bnd_flux_dn_dir
        1503             :
        1504     6602436 :       do i = 1, nday
        1505    17568900 :          cam_out%soll(idxday(i)) = sum(fsw%bnd_flux_dn_dir(i,nlay+1,1:9))      &
        1506    76131900 :                                    + 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10)
        1507             :
        1508     5856300 :          cam_out%sols(idxday(i)) = 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10)   &
        1509    35137800 :                                    + sum(fsw%bnd_flux_dn_dir(i,nlay+1,11:14))
        1510             :
        1511    11712600 :          cam_out%solld(idxday(i)) = sum(flux_dn_diffuse(i,nlay+1,1:9))         &
        1512    70275600 :                                     + 0.5_r8 * flux_dn_diffuse(i,nlay+1,10)
        1513             :
        1514     5856300 :          cam_out%solsd(idxday(i)) = 0.5_r8 * flux_dn_diffuse(i, nlay+1, 10)    &
        1515    41740236 :                                     + sum(flux_dn_diffuse(i,nlay+1,11:14))
        1516             :       end do
        1517             :
        1518     1492272 :    end subroutine set_sw_diags
        1519             :
        1520             :    !-------------------------------------------------------------------------------
        1521             :
        1522      746136 :    subroutine set_lw_diags()
        1523             :
        1524             :       ! Set CAM LW diagnostics
        1525             :       !----------------------------------------------------------------------------
        1526             :
        1527      746136 :       fnl = 0._r8
        1528      746136 :       fcnl = 0._r8
        1529             :
        1530             :       ! RTE-RRTMGP convention for net is (down - up) **CAM assumes (up - down) !!
        1531  1171867320 :       fnl(:ncol,ktopcam:)  = -1._r8 * flw%flux_net(    :, ktoprad:)
        1532  1171867320 :       fcnl(:ncol,ktopcam:) = -1._r8 * flwc%flux_net(   :, ktoprad:)
        1533             :
        1534  1171867320 :       rd%flux_lw_up(:ncol,ktopcam:)     = flw%flux_up( :, ktoprad:)
        1535  1171867320 :       rd%flux_lw_clr_up(:ncol,ktopcam:) = flwc%flux_up(:, ktoprad:)
        1536  1171867320 :       rd%flux_lw_dn(:ncol,ktopcam:)     = flw%flux_dn( :, ktoprad:)
        1537  1171867320 :       rd%flux_lw_clr_dn(:ncol,ktopcam:) = flwc%flux_dn(:, ktoprad:)
        1538             :
        1539      746136 :       call heating_rate('LW', ncol, fnl, qrl)
        1540      746136 :       call heating_rate('LW', ncol, fcnl, rd%qrlc)
        1541             :
        1542    12458736 :       flns(:ncol) = fnl(:ncol, pverp)
        1543    12458736 :       flnt(:ncol) = fnl(:ncol, ktopcam)
        1544             :
        1545    13204872 :       rd%flnsc(:ncol) = fcnl(:ncol, pverp)
        1546    13204872 :       rd%flntc(:ncol) = fcnl(:ncol, ktopcam)    ! net lw flux at top-of-model
        1547             :
        1548    12458736 :       cam_out%flwds(:ncol) = flw%flux_dn(:, nlay+1)
        1549    12458736 :       rd%fldsc(:ncol)      = flwc%flux_dn(:, nlay+1)
        1550             :
        1551    12458736 :       rd%flut(:ncol)  = flw%flux_up(:, ktoprad)
        1552    12458736 :       rd%flutc(:ncol) = flwc%flux_up(:, ktoprad)
        1553             :
        1554             :       ! Output fluxes at 200 mb
        1555      746136 :       call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl,  rd%fln200)
        1556      746136 :       call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
        1557      746136 :       if (hist_fld_active('FLNR')) then
        1558           0 :          do i = 1,ncol
        1559           0 :             call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
        1560             :          end do
        1561             :       end if
        1562             :
        1563      746136 :       if (spectralflux) then
        1564           0 :          lu  = 0._r8
        1565           0 :          ld  = 0._r8
        1566           0 :          lu(:ncol, ktopcam:, :)  = flw%bnd_flux_up(:, ktoprad:, :)
        1567           0 :          ld(:ncol, ktopcam:, :)  = flw%bnd_flux_dn(:, ktoprad:, :)
        1568             :       end if
        1569             :
        1570      746136 :    end subroutine set_lw_diags
        1571             :
        1572             :    !-------------------------------------------------------------------------------
        1573             :
        1574     2984544 :    subroutine heating_rate(type, ncol, flux_net, hrate)
        1575             :
        1576             :       ! Compute heating rate as a dry static energy tendency
        1577             :
        1578             :       ! arguments
        1579             :       character(2), intent(in)  :: type ! either LW or SW
        1580             :       integer,      intent(in)  :: ncol
        1581             :       real(r8),     intent(in)  :: flux_net(pcols,pverp)  ! W/m^2
        1582             :       real(r8),     intent(out) :: hrate(pcols,pver)      ! J/kg/s
        1583             :
        1584             :       ! local vars
        1585             :       integer :: k
        1586             :
        1587             :       ! Initialize for layers where RRTMGP is not providing fluxes.
        1588     2984544 :       hrate = 0.0_r8
        1589             :
        1590     1492272 :       select case (type)
        1591             :       case ('LW')
        1592             :
        1593   140273568 :          do k = ktopcam, pver
        1594             :             ! (flux divergence as bottom-MINUS-top) * g/dp
        1595   277562592 :             hrate(:ncol,k) = (flux_net(:ncol,k+1) - flux_net(:ncol,k)) * &
        1596  2596379760 :                           gravit * state%rpdel(:ncol,k)
        1597             :          end do
        1598             :
        1599             :       case ('SW')
        1600             :
        1601   141765840 :          do k = ktopcam, pver
        1602             :             ! top - bottom
        1603   277562592 :             hrate(:ncol,k) = (flux_net(:ncol,k) - flux_net(:ncol,k+1)) * &
        1604  2596379760 :                           gravit * state%rpdel(:ncol,k)
        1605             :          end do
        1606             :
        1607             :       end select
        1608             :
        1609     2984544 :    end subroutine heating_rate
        1610             :
        1611             :    !----------------------------------------------------------------------------
        1612             :    !            -- end contains statement of radiation_tend --
        1613             :    !----------------------------------------------------------------------------
        1614             : end subroutine radiation_tend
        1615             :
        1616             : !===============================================================================
        1617             :
        1618      746136 : subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
        1619             :
        1620             :    ! Dump shortwave radiation information to history buffer.
        1621             :
        1622             :    integer ,               intent(in) :: lchnk
        1623             :    integer,                intent(in) :: ncol
        1624             :    integer,                intent(in) :: icall
        1625             :    type(rad_out_t),        intent(in) :: rd
        1626             :    type(physics_buffer_desc), pointer :: pbuf(:)
        1627             :    type(cam_out_t),        intent(in) :: cam_out
        1628             :
        1629             :    ! local variables
        1630      746136 :    real(r8), pointer :: qrs(:,:)
        1631      746136 :    real(r8), pointer :: fsnt(:)
        1632      746136 :    real(r8), pointer :: fsns(:)
        1633      746136 :    real(r8), pointer :: fsds(:)
        1634             :
        1635             :    real(r8) :: ftem(pcols)
        1636             :    !----------------------------------------------------------------------------
        1637             :
        1638      746136 :    call pbuf_get_field(pbuf, qrs_idx,  qrs)
        1639      746136 :    call pbuf_get_field(pbuf, fsnt_idx, fsnt)
        1640      746136 :    call pbuf_get_field(pbuf, fsns_idx, fsns)
        1641      746136 :    call pbuf_get_field(pbuf, fsds_idx, fsds)
        1642             :
        1643      746136 :    call outfld('SOLIN'//diag(icall),    rd%solin,      pcols, lchnk)
        1644             :
        1645             :    ! QRS is output as temperature tendency.
        1646  1159408584 :    call outfld('QRS'//diag(icall),      qrs(:ncol,:)/cpair,     ncol, lchnk)
        1647  1159408584 :    call outfld('QRSC'//diag(icall),     rd%qrsc(:ncol,:)/cpair, ncol, lchnk)
        1648             :
        1649      746136 :    call outfld('FSNT'//diag(icall),    fsnt,               pcols, lchnk)
        1650      746136 :    call outfld('FSNTC'//diag(icall),   rd%fsntc,           pcols, lchnk)
        1651      746136 :    call outfld('FSNTOA'//diag(icall),  rd%fsntoa,          pcols, lchnk)
        1652      746136 :    call outfld('FSNTOAC'//diag(icall), rd%fsntoac,         pcols, lchnk)
        1653             :
        1654    12458736 :    ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol)
        1655      746136 :    call outfld('SWCF'//diag(icall),     ftem,          pcols, lchnk)
        1656             :
        1657      746136 :    call outfld('FSUTOA'//diag(icall),   rd%fsutoa,     pcols, lchnk)
        1658             :
        1659      746136 :    call outfld('FSNIRTOA'//diag(icall), rd%fsnirt,     pcols, lchnk)
        1660      746136 :    call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc,     pcols, lchnk)
        1661      746136 :    call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq,   pcols, lchnk)
        1662             :
        1663      746136 :    call outfld('FSN200'//diag(icall),   rd%fsn200,     pcols, lchnk)
        1664      746136 :    call outfld('FSN200C'//diag(icall),  rd%fsn200c,    pcols, lchnk)
        1665             :
        1666      746136 :    call outfld('FSNR'//diag(icall),     rd%fsnr,       pcols, lchnk)
        1667             :
        1668      746136 :    call outfld('SOLS'//diag(icall),     cam_out%sols,  pcols, lchnk)
        1669      746136 :    call outfld('SOLL'//diag(icall),     cam_out%soll,  pcols, lchnk)
        1670      746136 :    call outfld('SOLSD'//diag(icall),    cam_out%solsd, pcols, lchnk)
        1671      746136 :    call outfld('SOLLD'//diag(icall),    cam_out%solld, pcols, lchnk)
        1672             :
        1673      746136 :    call outfld('FSNS'//diag(icall),     fsns,          pcols, lchnk)
        1674      746136 :    call outfld('FSNSC'//diag(icall),    rd%fsnsc,      pcols, lchnk)
        1675             :
        1676      746136 :    call outfld('FSDS'//diag(icall),     fsds,          pcols, lchnk)
        1677      746136 :    call outfld('FSDSC'//diag(icall),    rd%fsdsc,      pcols, lchnk)
        1678             :
        1679      746136 :    call outfld('FUS'//diag(icall),  rd%flux_sw_up,     pcols, lchnk)
        1680      746136 :    call outfld('FUSC'//diag(icall), rd%flux_sw_clr_up, pcols, lchnk)
        1681      746136 :    call outfld('FDS'//diag(icall),  rd%flux_sw_dn,     pcols, lchnk)
        1682      746136 :    call outfld('FDSC'//diag(icall), rd%flux_sw_clr_dn, pcols, lchnk)
        1683             :
        1684      746136 : end subroutine radiation_output_sw
        1685             :
        1686             : !===============================================================================
        1687             :
        1688      746136 : subroutine radiation_output_cld(lchnk, rd)
        1689             :
        1690             :    ! Dump shortwave cloud optics information to history buffer.
        1691             :
        1692             :    integer ,        intent(in) :: lchnk
        1693             :    type(rad_out_t), intent(in) :: rd
        1694             :    !----------------------------------------------------------------------------
        1695             :
        1696      746136 :    call outfld('TOT_CLD_VISTAU',  rd%tot_cld_vistau,  pcols, lchnk)
        1697      746136 :    call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk)
        1698      746136 :    call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk)
        1699      746136 :    call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk)
        1700      746136 :    if (cldfsnow_idx > 0) then
        1701      746136 :       call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk)
        1702             :    endif
        1703      746136 :    if (cldfgrau_idx > 0 .and. graupel_in_rad) then
        1704           0 :       call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau , pcols, lchnk)
        1705             :    endif
        1706             :
        1707      746136 : end subroutine radiation_output_cld
        1708             :
        1709             : !===============================================================================
        1710             :
        1711      746136 : subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
        1712             :
        1713             :    ! Dump longwave radiation information to history buffer
        1714             :
        1715             :    integer,                intent(in) :: lchnk
        1716             :    integer,                intent(in) :: ncol
        1717             :    integer,                intent(in) :: icall  ! icall=0 for climate diagnostics
        1718             :    type(rad_out_t),        intent(in) :: rd
        1719             :    type(physics_buffer_desc), pointer :: pbuf(:)
        1720             :    type(cam_out_t),        intent(in) :: cam_out
        1721             :
        1722             :    ! local variables
        1723      746136 :    real(r8), pointer :: qrl(:,:)
        1724      746136 :    real(r8), pointer :: flnt(:)
        1725      746136 :    real(r8), pointer :: flns(:)
        1726             :
        1727             :    real(r8) :: ftem(pcols)
        1728             :    !----------------------------------------------------------------------------
        1729             :
        1730      746136 :    call pbuf_get_field(pbuf, qrl_idx,  qrl)
        1731      746136 :    call pbuf_get_field(pbuf, flnt_idx, flnt)
        1732      746136 :    call pbuf_get_field(pbuf, flns_idx, flns)
        1733             :
        1734  1159408584 :    call outfld('QRL'//diag(icall),     qrl(:ncol,:)/cpair,     ncol, lchnk)
        1735  1159408584 :    call outfld('QRLC'//diag(icall),    rd%qrlc(:ncol,:)/cpair, ncol, lchnk)
        1736             :
        1737      746136 :    call outfld('FLNT'//diag(icall),    flnt,          pcols, lchnk)
        1738      746136 :    call outfld('FLNTC'//diag(icall),   rd%flntc,      pcols, lchnk)
        1739             :
        1740      746136 :    call outfld('FLUT'//diag(icall),    rd%flut,       pcols, lchnk)
        1741      746136 :    call outfld('FLUTC'//diag(icall),   rd%flutc,      pcols, lchnk)
        1742             :
        1743    12458736 :    ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol)
        1744      746136 :    call outfld('LWCF'//diag(icall),    ftem,          pcols, lchnk)
        1745             :
        1746      746136 :    call outfld('FLN200'//diag(icall),  rd%fln200,     pcols, lchnk)
        1747      746136 :    call outfld('FLN200C'//diag(icall), rd%fln200c,    pcols, lchnk)
        1748             :
        1749      746136 :    call outfld('FLNR'//diag(icall),    rd%flnr,       pcols, lchnk)
        1750             :
        1751      746136 :    call outfld('FLNS'//diag(icall),    flns,          pcols, lchnk)
        1752      746136 :    call outfld('FLNSC'//diag(icall),   rd%flnsc,      pcols, lchnk)
        1753             :
        1754      746136 :    call outfld('FLDS'//diag(icall),    cam_out%flwds, pcols, lchnk)
        1755      746136 :    call outfld('FLDSC'//diag(icall),   rd%fldsc,      pcols, lchnk)
        1756             :
        1757      746136 :    call outfld('FDL'//diag(icall),  rd%flux_lw_dn,     pcols, lchnk)
        1758      746136 :    call outfld('FDLC'//diag(icall), rd%flux_lw_clr_dn, pcols, lchnk)
        1759      746136 :    call outfld('FUL'//diag(icall),  rd%flux_lw_up,     pcols, lchnk)
        1760      746136 :    call outfld('FULC'//diag(icall), rd%flux_lw_clr_up, pcols, lchnk)
        1761             :
        1762      746136 : end subroutine radiation_output_lw
        1763             :
        1764             : !===============================================================================
        1765             :
        1766        3072 : subroutine coefs_init(coefs_file, available_gases, kdist)
        1767             :
        1768             :    ! Read data from coefficients file.  Initialize the kdist object.
        1769             :    ! available_gases object provides the gas names that CAM provides.
        1770             :
        1771             :    ! arguments
        1772             :    character(len=*),                  intent(in)  :: coefs_file
        1773             :    class(ty_gas_concs),               intent(in)  :: available_gases
        1774             :    class(ty_gas_optics_rrtmgp),       intent(out) :: kdist
        1775             :
        1776             :    ! local variables
        1777             :    type(file_desc_t) :: fh    ! pio file handle
        1778             :    character(len=cl) :: locfn ! path to file on local storage
        1779             :
        1780             :    ! File dimensions
        1781             :    integer ::            &
        1782             :       absorber,          &
        1783             :       atmos_layer,       &
        1784             :       bnd,               &
        1785             :       pressure,          &
        1786             :       temperature,       &
        1787             :       absorber_ext,      &
        1788             :       pressure_interp,   &
        1789             :       mixing_fraction,   &
        1790             :       gpt,               &
        1791             :       temperature_Planck
        1792             :
        1793             :    integer :: i
        1794             :    integer :: did, vid
        1795             :    integer :: ierr, istat
        1796             :
        1797        3072 :    character(32), dimension(:),  allocatable :: gas_names
        1798        3072 :    integer,  dimension(:,:,:),   allocatable :: key_species
        1799        3072 :    integer,  dimension(:,:),     allocatable :: band2gpt
        1800        3072 :    real(r8), dimension(:,:),     allocatable :: band_lims_wavenum
        1801        3072 :    real(r8), dimension(:),       allocatable :: press_ref, temp_ref
        1802             :    real(r8)                                  :: press_ref_trop, temp_ref_t, temp_ref_p
        1803        3072 :    real(r8), dimension(:,:,:),   allocatable :: vmr_ref
        1804        3072 :    real(r8), dimension(:,:,:,:), allocatable :: kmajor
        1805        3072 :    real(r8), dimension(:,:,:),   allocatable :: kminor_lower, kminor_upper
        1806        3072 :    real(r8), dimension(:,:),     allocatable :: totplnk
        1807        3072 :    real(r8), dimension(:,:,:,:), allocatable :: planck_frac
        1808        3072 :    real(r8), dimension(:),       allocatable :: solar_src_quiet, solar_src_facular, solar_src_sunspot
        1809             :    real(r8)                                  :: tsi_default
        1810        3072 :    real(r8), dimension(:,:,:),   allocatable :: rayl_lower, rayl_upper
        1811        3072 :    character(len=32), dimension(:),  allocatable :: gas_minor,         &
        1812        3072 :                                                     identifier_minor,  &
        1813        3072 :                                                     minor_gases_lower, &
        1814        3072 :                                                     minor_gases_upper, &
        1815        3072 :                                                     scaling_gas_lower, &
        1816        3072 :                                                     scaling_gas_upper
        1817        3072 :    integer, dimension(:,:),          allocatable :: minor_limits_gpt_lower, &
        1818        3072 :                                                     minor_limits_gpt_upper
        1819             :    ! Send these to RRTMGP as logicals,
        1820             :    ! but they have to be read from the netCDF as integers
        1821        3072 :    logical, dimension(:),            allocatable :: minor_scales_with_density_lower, &
        1822        3072 :                                                     minor_scales_with_density_upper
        1823        3072 :    logical, dimension(:),            allocatable :: scale_by_complement_lower, &
        1824        3072 :                                                     scale_by_complement_upper
        1825        3072 :    integer, dimension(:), allocatable :: int2log   ! use this to convert integer-to-logical.
        1826        3072 :    integer, dimension(:),            allocatable :: kminor_start_lower, kminor_start_upper
        1827        3072 :    real(r8), dimension(:,:),         allocatable :: optimal_angle_fit
        1828             :    real(r8)                                      :: mg_default, sb_default
        1829             :
        1830             :    integer :: pairs, &
        1831             :               minorabsorbers, &
        1832             :               minor_absorber_intervals_lower, &
        1833             :               minor_absorber_intervals_upper, &
        1834             :               contributors_lower, &
        1835             :               contributors_upper, &
        1836             :               fit_coeffs
        1837             :
        1838             :    character(len=128) :: error_msg
        1839             :    character(len=*), parameter :: sub = 'coefs_init'
        1840             :    !----------------------------------------------------------------------------
        1841             :
        1842             :    ! Open file
        1843        3072 :    call getfil(coefs_file, locfn, 0)
        1844        3072 :    call cam_pio_openfile(fh, locfn, PIO_NOWRITE)
        1845             :
        1846        3072 :    call pio_seterrorhandling(fh, PIO_BCAST_ERROR)
        1847             :
        1848             :    ! Get dimensions
        1849             :
        1850        3072 :    ierr = pio_inq_dimid(fh, 'absorber', did)
        1851        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': absorber not found')
        1852        3072 :    ierr = pio_inq_dimlen(fh, did, absorber)
        1853             :
        1854        3072 :    ierr = pio_inq_dimid(fh, 'atmos_layer', did)
        1855        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': atmos_layer not found')
        1856        3072 :    ierr = pio_inq_dimlen(fh, did, atmos_layer)
        1857             :
        1858        3072 :    ierr = pio_inq_dimid(fh, 'bnd', did)
        1859        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': bnd not found')
        1860        3072 :    ierr = pio_inq_dimlen(fh, did, bnd)
        1861             :
        1862        3072 :    ierr = pio_inq_dimid(fh, 'pressure', did)
        1863        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': pressure not found')
        1864        3072 :    ierr = pio_inq_dimlen(fh, did, pressure)
        1865             :
        1866        3072 :    ierr = pio_inq_dimid(fh, 'temperature', did)
        1867        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': temperature not found')
        1868        3072 :    ierr = pio_inq_dimlen(fh, did, temperature)
        1869             :
        1870        3072 :    ierr = pio_inq_dimid(fh, 'absorber_ext', did)
        1871        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': absorber_ext not found')
        1872        3072 :    ierr = pio_inq_dimlen(fh, did, absorber_ext)
        1873             :
        1874        3072 :    ierr = pio_inq_dimid(fh, 'pressure_interp', did)
        1875        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': pressure_interp not found')
        1876        3072 :    ierr = pio_inq_dimlen(fh, did, pressure_interp)
        1877             :
        1878        3072 :    ierr = pio_inq_dimid(fh, 'mixing_fraction', did)
        1879        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': mixing_fraction not found')
        1880        3072 :    ierr = pio_inq_dimlen(fh, did, mixing_fraction)
        1881             :
        1882        3072 :    ierr = pio_inq_dimid(fh, 'gpt', did)
        1883        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': gpt not found')
        1884        3072 :    ierr = pio_inq_dimlen(fh, did, gpt)
        1885             :
        1886        3072 :    temperature_Planck = 0
        1887        3072 :    ierr = pio_inq_dimid(fh, 'temperature_Planck', did)
        1888        3072 :    if (ierr == PIO_NOERR) then
        1889        1536 :       ierr = pio_inq_dimlen(fh, did, temperature_Planck)
        1890             :    end if
        1891        3072 :    ierr = pio_inq_dimid(fh, 'pair', did)
        1892        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': pair not found')
        1893        3072 :    ierr = pio_inq_dimlen(fh, did, pairs)
        1894        3072 :    ierr = pio_inq_dimid(fh, 'minor_absorber', did)
        1895        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber not found')
        1896        3072 :    ierr = pio_inq_dimlen(fh, did, minorabsorbers)
        1897        3072 :    ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_lower', did)
        1898        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_lower not found')
        1899        3072 :    ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_lower)
        1900        3072 :    ierr = pio_inq_dimid(fh, 'minor_absorber_intervals_upper', did)
        1901        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_absorber_intervals_upper not found')
        1902        3072 :    ierr = pio_inq_dimlen(fh, did, minor_absorber_intervals_upper)
        1903        3072 :    ierr = pio_inq_dimid(fh, 'contributors_lower', did)
        1904        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': contributors_lower not found')
        1905        3072 :    ierr = pio_inq_dimlen(fh, did, contributors_lower)
        1906        3072 :    ierr = pio_inq_dimid(fh, 'contributors_upper', did)
        1907        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': contributors_upper not found')
        1908        3072 :    ierr = pio_inq_dimlen(fh, did, contributors_upper)
        1909             :
        1910        3072 :    ierr = pio_inq_dimid(fh, 'fit_coeffs', did)
        1911        3072 :    if (ierr == PIO_NOERR) then
        1912        1536 :       ierr = pio_inq_dimlen(fh, did, fit_coeffs)
        1913             :    end if
        1914             :
        1915             :    ! Get variables
        1916             :
        1917             :    ! names of absorbing gases
        1918        9216 :    allocate(gas_names(absorber), stat=istat)
        1919        3072 :    call handle_allocate_error(istat, sub, 'gas_names')
        1920        3072 :    ierr = pio_inq_varid(fh, 'gas_names', vid)
        1921        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': gas_names not found')
        1922        3072 :    ierr = pio_get_var(fh, vid, gas_names)
        1923        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_names')
        1924             :
        1925             :    ! key species pair for each band
        1926       12288 :    allocate(key_species(2,atmos_layer,bnd), stat=istat)
        1927        3072 :    call handle_allocate_error(istat, sub, 'key_species')
        1928        3072 :    ierr = pio_inq_varid(fh, 'key_species', vid)
        1929        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': key_species not found')
        1930        3072 :    ierr = pio_get_var(fh, vid, key_species)
        1931        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading key_species')
        1932             :
        1933             :    ! beginning and ending gpoint for each band
        1934        9216 :    allocate(band2gpt(2,bnd), stat=istat)
        1935        3072 :    call handle_allocate_error(istat, sub, 'band2gpt')
        1936        3072 :    ierr = pio_inq_varid(fh, 'bnd_limits_gpt', vid)
        1937        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_gpt not found')
        1938        3072 :    ierr = pio_get_var(fh, vid, band2gpt)
        1939        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_gpt')
        1940             :
        1941             :    ! beginning and ending wavenumber for each band
        1942        9216 :    allocate(band_lims_wavenum(2,bnd), stat=istat)
        1943        3072 :    call handle_allocate_error(istat, sub, 'band_lims_wavenum')
        1944        3072 :    ierr = pio_inq_varid(fh, 'bnd_limits_wavenumber', vid)
        1945        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': bnd_limits_wavenumber not found')
        1946        3072 :    ierr = pio_get_var(fh, vid, band_lims_wavenum)
        1947        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading bnd_limits_wavenumber')
        1948             :
        1949             :    ! pressures [hPa] for reference atmosphere; press_ref(# reference layers)
        1950        9216 :    allocate(press_ref(pressure), stat=istat)
        1951        3072 :    call handle_allocate_error(istat, sub, 'press_ref')
        1952        3072 :    ierr = pio_inq_varid(fh, 'press_ref', vid)
        1953        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': press_ref not found')
        1954        3072 :    ierr = pio_get_var(fh, vid, press_ref)
        1955        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref')
        1956             :
        1957             :    ! reference pressure separating the lower and upper atmosphere
        1958        3072 :    ierr = pio_inq_varid(fh, 'press_ref_trop', vid)
        1959        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': press_ref_trop not found')
        1960        3072 :    ierr = pio_get_var(fh, vid, press_ref_trop)
        1961        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading press_ref_trop')
        1962             :
        1963             :    ! temperatures [K] for reference atmosphere; temp_ref(# reference layers)
        1964        9216 :    allocate(temp_ref(temperature), stat=istat)
        1965        3072 :    call handle_allocate_error(istat, sub, 'temp_ref')
        1966        3072 :    ierr = pio_inq_varid(fh, 'temp_ref', vid)
        1967        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': temp_ref not found')
        1968        3072 :    ierr = pio_get_var(fh, vid, temp_ref)
        1969        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading temp_ref')
        1970             :
        1971             :    ! standard spectroscopic reference temperature [K]
        1972        3072 :    ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_T', vid)
        1973        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_T not found')
        1974        3072 :    ierr = pio_get_var(fh, vid, temp_ref_t)
        1975        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_T')
        1976             :
        1977             :    ! standard spectroscopic reference pressure [hPa]
        1978        3072 :    ierr = pio_inq_varid(fh, 'absorption_coefficient_ref_P', vid)
        1979        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': absorption_coefficient_ref_P not found')
        1980        3072 :    ierr = pio_get_var(fh, vid, temp_ref_p)
        1981        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading absorption_coefficient_ref_P')
        1982             :
        1983             :    ! volume mixing ratios for reference atmosphere
        1984       15360 :    allocate(vmr_ref(atmos_layer, absorber_ext, temperature), stat=istat)
        1985        3072 :    call handle_allocate_error(istat, sub, 'vmr_ref')
        1986        3072 :    ierr = pio_inq_varid(fh, 'vmr_ref', vid)
        1987        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': vmr_ref not found')
        1988        3072 :    ierr = pio_get_var(fh, vid, vmr_ref)
        1989        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading vmr_ref')
        1990             :
        1991             :    ! absorption coefficients due to major absorbing gases
        1992       18432 :    allocate(kmajor(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
        1993        3072 :    call handle_allocate_error(istat, sub, 'kmajor')
        1994        3072 :    ierr = pio_inq_varid(fh, 'kmajor', vid)
        1995        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': kmajor not found')
        1996        3072 :    ierr = pio_get_var(fh, vid, kmajor)
        1997        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading kmajor')
        1998             :
        1999             :    ! absorption coefficients due to minor absorbing gases in lower part of atmosphere
        2000       15360 :    allocate(kminor_lower(contributors_lower, mixing_fraction, temperature), stat=istat)
        2001        3072 :    call handle_allocate_error(istat, sub, 'kminor_lower')
        2002        3072 :    ierr = pio_inq_varid(fh, 'kminor_lower', vid)
        2003        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': kminor_lower not found')
        2004        3072 :    ierr = pio_get_var(fh, vid, kminor_lower)
        2005        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_lower')
        2006             :
        2007             :    ! absorption coefficients due to minor absorbing gases in upper part of atmosphere
        2008       15360 :    allocate(kminor_upper(contributors_upper, mixing_fraction, temperature), stat=istat)
        2009        3072 :    call handle_allocate_error(istat, sub, 'kminor_upper')
        2010        3072 :    ierr = pio_inq_varid(fh, 'kminor_upper', vid)
        2011        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': kminor_upper not found')
        2012        3072 :    ierr = pio_get_var(fh, vid, kminor_upper)
        2013        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_upper')
        2014             :
        2015             :    ! integrated Planck function by band
        2016        3072 :    ierr = pio_inq_varid(fh, 'totplnk', vid)
        2017        3072 :    if (ierr == PIO_NOERR) then
        2018        6144 :       allocate(totplnk(temperature_Planck,bnd), stat=istat)
        2019        1536 :       call handle_allocate_error(istat, sub, 'totplnk')
        2020        1536 :       ierr = pio_get_var(fh, vid, totplnk)
        2021        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading totplnk')
        2022             :    end if
        2023             :
        2024             :    ! Planck fractions
        2025        3072 :    ierr = pio_inq_varid(fh, 'plank_fraction', vid)
        2026        3072 :    if (ierr == PIO_NOERR) then
        2027        9216 :       allocate(planck_frac(gpt,mixing_fraction,pressure_interp,temperature), stat=istat)
        2028        1536 :       call handle_allocate_error(istat, sub, 'planck_frac')
        2029        1536 :       ierr = pio_get_var(fh, vid, planck_frac)
        2030        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading plank_fraction')
        2031             :    end if
        2032             :
        2033        3072 :    ierr = pio_inq_varid(fh, 'optimal_angle_fit', vid)
        2034        3072 :    if (ierr == PIO_NOERR) then
        2035        6144 :       allocate(optimal_angle_fit(fit_coeffs, bnd), stat=istat)
        2036        1536 :       call handle_allocate_error(istat, sub, 'optiman_angle_fit')
        2037        1536 :       ierr = pio_get_var(fh, vid, optimal_angle_fit)
        2038        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading optimal_angle_fit')
        2039             :    end if
        2040             :
        2041        3072 :    ierr = pio_inq_varid(fh, 'solar_source_quiet', vid)
        2042        3072 :    if (ierr == PIO_NOERR) then
        2043        4608 :       allocate(solar_src_quiet(gpt), stat=istat)
        2044        1536 :       call handle_allocate_error(istat, sub, 'solar_src_quiet')
        2045        1536 :       ierr = pio_get_var(fh, vid, solar_src_quiet)
        2046        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_quiet')
        2047             :    end if
        2048             :
        2049        3072 :    ierr = pio_inq_varid(fh, 'solar_source_facular', vid)
        2050        3072 :    if (ierr == PIO_NOERR) then
        2051        4608 :       allocate(solar_src_facular(gpt), stat=istat)
        2052        1536 :       call handle_allocate_error(istat, sub, 'solar_src_facular')
        2053        1536 :       ierr = pio_get_var(fh, vid, solar_src_facular)
        2054        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_facular')
        2055             :    end if
        2056             :
        2057        3072 :    ierr = pio_inq_varid(fh, 'solar_source_sunspot', vid)
        2058        3072 :    if (ierr == PIO_NOERR) then
        2059        4608 :       allocate(solar_src_sunspot(gpt), stat=istat)
        2060        1536 :       call handle_allocate_error(istat, sub, 'solar_src_sunspot')
        2061        1536 :       ierr = pio_get_var(fh, vid, solar_src_sunspot)
        2062        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_sunspot')
        2063             :    end if
        2064             :
        2065        3072 :    ierr = pio_inq_varid(fh, 'tsi_default', vid)
        2066        3072 :    if (ierr == PIO_NOERR) then
        2067        1536 :       ierr = pio_get_var(fh, vid, tsi_default)
        2068        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading tsi_default')
        2069             :    end if
        2070             :
        2071        3072 :    ierr = pio_inq_varid(fh, 'mg_default', vid)
        2072        3072 :    if (ierr == PIO_NOERR) then
        2073        1536 :       ierr = pio_get_var(fh, vid, mg_default)
        2074        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading mg_default')
        2075             :    end if
        2076             :
        2077        3072 :    ierr = pio_inq_varid(fh, 'sb_default', vid)
        2078        3072 :    if (ierr == PIO_NOERR) then
        2079        1536 :       ierr = pio_get_var(fh, vid, sb_default)
        2080        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading sb_default')
        2081             :    end if
        2082             :
        2083             :    ! rayleigh scattering contribution in lower part of atmosphere
        2084        3072 :    ierr = pio_inq_varid(fh, 'rayl_lower', vid)
        2085        3072 :    if (ierr == PIO_NOERR) then
        2086        7680 :       allocate(rayl_lower(gpt,mixing_fraction,temperature), stat=istat)
        2087        1536 :       call handle_allocate_error(istat, sub, 'rayl_lower')
        2088        1536 :       ierr = pio_get_var(fh, vid, rayl_lower)
        2089        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_lower')
        2090             :    end if
        2091             :
        2092             :    ! rayleigh scattering contribution in upper part of atmosphere
        2093        3072 :    ierr = pio_inq_varid(fh, 'rayl_upper', vid)
        2094        3072 :    if (ierr == PIO_NOERR) then
        2095        7680 :       allocate(rayl_upper(gpt,mixing_fraction,temperature), stat=istat)
        2096        1536 :       call handle_allocate_error(istat, sub, 'rayl_upper')
        2097        1536 :       ierr = pio_get_var(fh, vid, rayl_upper)
        2098        1536 :       if (ierr /= PIO_NOERR) call endrun(sub//': error reading rayl_upper')
        2099             :    end if
        2100             :
        2101        9216 :    allocate(gas_minor(minorabsorbers), stat=istat)
        2102        3072 :    call handle_allocate_error(istat, sub, 'gas_minor')
        2103        3072 :    ierr = pio_inq_varid(fh, 'gas_minor', vid)
        2104        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': gas_minor not found')
        2105        3072 :    ierr = pio_get_var(fh, vid, gas_minor)
        2106        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading gas_minor')
        2107             :
        2108        9216 :    allocate(identifier_minor(minorabsorbers), stat=istat)
        2109        3072 :    call handle_allocate_error(istat, sub, 'identifier_minor')
        2110        3072 :    ierr = pio_inq_varid(fh, 'identifier_minor', vid)
        2111        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': identifier_minor not found')
        2112        3072 :    ierr = pio_get_var(fh, vid, identifier_minor)
        2113        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading identifier_minor')
        2114             :
        2115        9216 :    allocate(minor_gases_lower(minor_absorber_intervals_lower), stat=istat)
        2116        3072 :    call handle_allocate_error(istat, sub, 'minor_gases_lower')
        2117        3072 :    ierr = pio_inq_varid(fh, 'minor_gases_lower', vid)
        2118        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_lower not found')
        2119        3072 :    ierr = pio_get_var(fh, vid, minor_gases_lower)
        2120        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_lower')
        2121             :
        2122        9216 :    allocate(minor_gases_upper(minor_absorber_intervals_upper), stat=istat)
        2123        3072 :    call handle_allocate_error(istat, sub, 'minor_gases_upper')
        2124        3072 :    ierr = pio_inq_varid(fh, 'minor_gases_upper', vid)
        2125        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_gases_upper not found')
        2126        3072 :    ierr = pio_get_var(fh, vid, minor_gases_upper)
        2127        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_gases_upper')
        2128             :
        2129       12288 :    allocate(minor_limits_gpt_lower(pairs,minor_absorber_intervals_lower), stat=istat)
        2130        3072 :    call handle_allocate_error(istat, sub, 'minor_limits_gpt_lower')
        2131        3072 :    ierr = pio_inq_varid(fh, 'minor_limits_gpt_lower', vid)
        2132        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_lower not found')
        2133        3072 :    ierr = pio_get_var(fh, vid, minor_limits_gpt_lower)
        2134        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_lower')
        2135             :
        2136       12288 :    allocate(minor_limits_gpt_upper(pairs,minor_absorber_intervals_upper), stat=istat)
        2137        3072 :    call handle_allocate_error(istat, sub, 'minor_limits_gpt_upper')
        2138        3072 :    ierr = pio_inq_varid(fh, 'minor_limits_gpt_upper', vid)
        2139        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_limits_gpt_upper not found')
        2140        3072 :    ierr = pio_get_var(fh, vid, minor_limits_gpt_upper)
        2141        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_limits_gpt_upper')
        2142             :
        2143             :    ! Read as integer and convert to logical
        2144        9216 :    allocate(int2log(minor_absorber_intervals_lower), stat=istat)
        2145        3072 :    call handle_allocate_error(istat, sub, 'int2log for lower')
        2146             :
        2147        9216 :    allocate(minor_scales_with_density_lower(minor_absorber_intervals_lower), stat=istat)
        2148        3072 :    call handle_allocate_error(istat, sub, 'minor_scales_with_density_lower')
        2149        3072 :    ierr = pio_inq_varid(fh, 'minor_scales_with_density_lower', vid)
        2150        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_lower not found')
        2151        3072 :    ierr = pio_get_var(fh, vid, int2log)
        2152        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_lower')
        2153      147456 :    do i = 1,minor_absorber_intervals_lower
        2154      147456 :       if (int2log(i) .eq. 0) then
        2155       53760 :          minor_scales_with_density_lower(i) = .false.
        2156             :       else
        2157       90624 :          minor_scales_with_density_lower(i) = .true.
        2158             :       end if
        2159             :    end do
        2160             :
        2161        9216 :    allocate(scale_by_complement_lower(minor_absorber_intervals_lower), stat=istat)
        2162        3072 :    call handle_allocate_error(istat, sub, 'scale_by_complement_lower')
        2163        3072 :    ierr = pio_inq_varid(fh, 'scale_by_complement_lower', vid)
        2164        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_lower not found')
        2165        3072 :    ierr = pio_get_var(fh, vid, int2log)
        2166        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_lower')
        2167      147456 :    do i = 1,minor_absorber_intervals_lower
        2168      147456 :       if (int2log(i) .eq. 0) then
        2169      104448 :          scale_by_complement_lower(i) = .false.
        2170             :       else
        2171       39936 :          scale_by_complement_lower(i) = .true.
        2172             :       end if
        2173             :    end do
        2174             :
        2175        3072 :    deallocate(int2log)
        2176             :
        2177             :    ! Read as integer and convert to logical
        2178        9216 :    allocate(int2log(minor_absorber_intervals_upper), stat=istat)
        2179        3072 :    call handle_allocate_error(istat, sub, 'int2log for upper')
        2180             :
        2181        9216 :    allocate(minor_scales_with_density_upper(minor_absorber_intervals_upper), stat=istat)
        2182        3072 :    call handle_allocate_error(istat, sub, 'minor_scales_with_density_upper')
        2183        3072 :    ierr = pio_inq_varid(fh, 'minor_scales_with_density_upper', vid)
        2184        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': minor_scales_with_density_upper not found')
        2185        3072 :    ierr = pio_get_var(fh, vid, int2log)
        2186        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading minor_scales_with_density_upper')
        2187       92160 :    do i = 1,minor_absorber_intervals_upper
        2188       92160 :       if (int2log(i) .eq. 0) then
        2189       46080 :          minor_scales_with_density_upper(i) = .false.
        2190             :       else
        2191       43008 :          minor_scales_with_density_upper(i) = .true.
        2192             :       end if
        2193             :    end do
        2194             :
        2195        9216 :    allocate(scale_by_complement_upper(minor_absorber_intervals_upper), stat=istat)
        2196        3072 :    call handle_allocate_error(istat, sub, 'scale_by_complement_upper')
        2197        3072 :    ierr = pio_inq_varid(fh, 'scale_by_complement_upper', vid)
        2198        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': scale_by_complement_upper not found')
        2199        3072 :    ierr = pio_get_var(fh, vid, int2log)
        2200        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading scale_by_complement_upper')
        2201       92160 :    do i = 1,minor_absorber_intervals_upper
        2202       92160 :       if (int2log(i) .eq. 0) then
        2203       70656 :          scale_by_complement_upper(i) = .false.
        2204             :       else
        2205       18432 :          scale_by_complement_upper(i) = .true.
        2206             :       end if
        2207             :    end do
        2208             :
        2209        3072 :    deallocate(int2log)
        2210             :
        2211        9216 :    allocate(scaling_gas_lower(minor_absorber_intervals_lower), stat=istat)
        2212        3072 :    call handle_allocate_error(istat, sub, 'scaling_gas_lower')
        2213        3072 :    ierr = pio_inq_varid(fh, 'scaling_gas_lower', vid)
        2214        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_lower not found')
        2215        3072 :    ierr = pio_get_var(fh, vid, scaling_gas_lower)
        2216        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_lower')
        2217             :
        2218        9216 :    allocate(scaling_gas_upper(minor_absorber_intervals_upper), stat=istat)
        2219        3072 :    call handle_allocate_error(istat, sub, 'scaling_gas_upper')
        2220        3072 :    ierr = pio_inq_varid(fh, 'scaling_gas_upper', vid)
        2221        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': scaling_gas_upper not found')
        2222        3072 :    ierr = pio_get_var(fh, vid, scaling_gas_upper)
        2223        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading scaling_gas_upper')
        2224             :
        2225        9216 :    allocate(kminor_start_lower(minor_absorber_intervals_lower), stat=istat)
        2226        3072 :    call handle_allocate_error(istat, sub, 'kminor_start_lower')
        2227        3072 :    ierr = pio_inq_varid(fh, 'kminor_start_lower', vid)
        2228        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_lower not found')
        2229        3072 :    ierr = pio_get_var(fh, vid, kminor_start_lower)
        2230        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_lower')
        2231             :
        2232        9216 :    allocate(kminor_start_upper(minor_absorber_intervals_upper), stat=istat)
        2233        3072 :    call handle_allocate_error(istat, sub, 'kminor_start_upper')
        2234        3072 :    ierr = pio_inq_varid(fh, 'kminor_start_upper', vid)
        2235        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': kminor_start_upper not found')
        2236        3072 :    ierr = pio_get_var(fh, vid, kminor_start_upper)
        2237        3072 :    if (ierr /= PIO_NOERR) call endrun(sub//': error reading kminor_start_upper')
        2238             :
        2239             :    ! Close file
        2240        3072 :    call pio_closefile(fh)
        2241             :
        2242             :    ! Initialize the gas optics object with data. The calls are slightly different depending
        2243             :    ! on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave)
        2244             :
        2245        3072 :    if (allocated(totplnk) .and. allocated(planck_frac)) then
        2246             :       error_msg = kdist%load( &
        2247             :          available_gases, gas_names, key_species,              &
        2248             :          band2gpt, band_lims_wavenum,                          &
        2249             :          press_ref, press_ref_trop, temp_ref,                  &
        2250             :          temp_ref_p, temp_ref_t, vmr_ref,                      &
        2251             :          kmajor, kminor_lower, kminor_upper,                   &
        2252             :          gas_minor, identifier_minor,                          &
        2253             :          minor_gases_lower, minor_gases_upper,                 &
        2254             :          minor_limits_gpt_lower, minor_limits_gpt_upper,       &
        2255             :          minor_scales_with_density_lower,                      &
        2256             :          minor_scales_with_density_upper,                      &
        2257             :          scaling_gas_lower, scaling_gas_upper,                 &
        2258             :          scale_by_complement_lower, scale_by_complement_upper, &
        2259             :          kminor_start_lower, kminor_start_upper,               &
        2260             :          totplnk, planck_frac, rayl_lower, rayl_upper,         &
        2261        1536 :          optimal_angle_fit)
        2262        1536 :    else if (allocated(solar_src_quiet)) then
        2263             :       error_msg = kdist%load( &
        2264             :          available_gases, gas_names, key_species,               &
        2265             :          band2gpt, band_lims_wavenum,                           &
        2266             :          press_ref, press_ref_trop, temp_ref,                   &
        2267             :          temp_ref_p, temp_ref_t, vmr_ref,                       &
        2268             :          kmajor, kminor_lower, kminor_upper,                    &
        2269             :          gas_minor, identifier_minor,                           &
        2270             :          minor_gases_lower, minor_gases_upper,                  &
        2271             :          minor_limits_gpt_lower, minor_limits_gpt_upper,        &
        2272             :          minor_scales_with_density_lower,                       &
        2273             :          minor_scales_with_density_upper,                       &
        2274             :          scaling_gas_lower, scaling_gas_upper,                  &
        2275             :          scale_by_complement_lower,                             &
        2276             :          scale_by_complement_upper,                             &
        2277             :          kminor_start_lower,                                    &
        2278             :          kminor_start_upper,                                    &
        2279             :          solar_src_quiet, solar_src_facular, solar_src_sunspot, &
        2280             :          tsi_default, mg_default, sb_default,                   &
        2281        1536 :          rayl_lower, rayl_upper)
        2282             :    else
        2283           0 :       error_msg = 'must supply either totplnk and planck_frac, or solar_src_[*]'
        2284             :    end if
        2285             :
        2286        3072 :    call stop_on_err(error_msg, sub, 'kdist%load')
        2287             :
        2288           0 :    deallocate( &
        2289           0 :       gas_names, key_species,               &
        2290           0 :       band2gpt, band_lims_wavenum,          &
        2291           0 :       press_ref, temp_ref, vmr_ref,         &
        2292           0 :       kmajor, kminor_lower, kminor_upper,   &
        2293           0 :       gas_minor, identifier_minor,          &
        2294           0 :       minor_gases_lower, minor_gases_upper, &
        2295           0 :       minor_limits_gpt_lower,               &
        2296           0 :       minor_limits_gpt_upper,               &
        2297           0 :       minor_scales_with_density_lower,      &
        2298           0 :       minor_scales_with_density_upper,      &
        2299           0 :       scale_by_complement_lower,            &
        2300           0 :       scale_by_complement_upper,            &
        2301           0 :       scaling_gas_lower, scaling_gas_upper, &
        2302        3072 :       kminor_start_lower, kminor_start_upper)
        2303             :
        2304        3072 :    if (allocated(totplnk))           deallocate(totplnk)
        2305        3072 :    if (allocated(planck_frac))       deallocate(planck_frac)
        2306        3072 :    if (allocated(optimal_angle_fit)) deallocate(optimal_angle_fit)
        2307        3072 :    if (allocated(solar_src_quiet))   deallocate(solar_src_quiet)
        2308        3072 :    if (allocated(solar_src_facular)) deallocate(solar_src_facular)
        2309        3072 :    if (allocated(solar_src_sunspot)) deallocate(solar_src_sunspot)
        2310        3072 :    if (allocated(rayl_lower))        deallocate(rayl_lower)
        2311        3072 :    if (allocated(rayl_upper))        deallocate(rayl_upper)
        2312             :
        2313        6144 : end subroutine coefs_init
        2314             :
        2315             : !=========================================================================================
        2316             :
        2317     5969088 : subroutine initialize_rrtmgp_fluxes(ncol, nlevels, nbands, fluxes, do_direct)
        2318             :
        2319             :    ! Allocate flux arrays and set values to zero.
        2320             :
        2321             :    ! Arguments
        2322             :    integer,                    intent(in)    :: ncol, nlevels, nbands
        2323             :    class(ty_fluxes_broadband), intent(inout) :: fluxes
        2324             :    logical, optional,          intent(in)    :: do_direct
        2325             :
        2326             :    ! Local variables
        2327             :    logical :: do_direct_local
        2328             :    integer :: istat
        2329             :    character(len=*), parameter :: sub = 'initialize_rrtmgp_fluxes'
        2330             :    !----------------------------------------------------------------------------
        2331             :
        2332     5969088 :    if (present(do_direct)) then
        2333             :       do_direct_local = .true.
        2334             :    else
        2335     2984544 :       do_direct_local = .false.
        2336             :    end if
        2337             :
        2338             :    ! Broadband fluxes
        2339    22448964 :    allocate(fluxes%flux_up(ncol, nlevels), stat=istat)
        2340     5969088 :    call handle_allocate_error(istat, sub, 'fluxes%flux_up')
        2341    16479876 :    allocate(fluxes%flux_dn(ncol, nlevels), stat=istat)
        2342     5969088 :    call handle_allocate_error(istat, sub, 'fluxes%flux_dn')
        2343    16479876 :    allocate(fluxes%flux_net(ncol, nlevels), stat=istat)
        2344     5969088 :    call handle_allocate_error(istat, sub, 'fluxes%flux_net')
        2345     5969088 :    if (do_direct_local) then
        2346     7526244 :       allocate(fluxes%flux_dn_dir(ncol, nlevels), stat=istat)
        2347     2984544 :       call handle_allocate_error(istat, sub, 'fluxes%flux_dn_dir')
        2348             :    end if
        2349             :
        2350             :    select type (fluxes)
        2351             :    type is (ty_fluxes_byband)
        2352             :       ! Fluxes by band always needed for SW.  Only allocate for LW
        2353             :       ! when spectralflux is true.
        2354     2984544 :       if (nbands == nswbands .or. spectralflux) then
        2355     6747666 :          allocate(fluxes%bnd_flux_up(ncol, nlevels, nbands), stat=istat)
        2356     1492272 :          call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_up')
        2357     5255394 :          allocate(fluxes%bnd_flux_dn(ncol, nlevels, nbands), stat=istat)
        2358     1492272 :          call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn')
        2359     5255394 :          allocate(fluxes%bnd_flux_net(ncol, nlevels, nbands), stat=istat)
        2360     1492272 :          call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_net')
        2361     1492272 :          if (do_direct_local) then
        2362     5255394 :             allocate(fluxes%bnd_flux_dn_dir(ncol, nlevels, nbands), stat=istat)
        2363     1492272 :             call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn_dir')
        2364             :          end if
        2365             :       end if
        2366             :    end select
        2367             :
        2368             :    ! Initialize
        2369     5969088 :    call reset_fluxes(fluxes)
        2370             :
        2371     5969088 : end subroutine initialize_rrtmgp_fluxes
        2372             :
        2373             : !=========================================================================================
        2374             :
        2375     5969088 : subroutine reset_fluxes(fluxes)
        2376             :
        2377             :    ! Reset flux arrays to zero.
        2378             :
        2379             :    class(ty_fluxes_broadband), intent(inout) :: fluxes
        2380             :    !----------------------------------------------------------------------------
        2381             :
        2382             :    ! Reset broadband fluxes
        2383  7249214448 :    fluxes%flux_up(:,:) = 0._r8
        2384  7249214448 :    fluxes%flux_dn(:,:) = 0._r8
        2385  7249214448 :    fluxes%flux_net(:,:) = 0._r8
        2386  2514894768 :    if (associated(fluxes%flux_dn_dir)) fluxes%flux_dn_dir(:,:) = 0._r8
        2387             :
        2388             :    select type (fluxes)
        2389             :    type is (ty_fluxes_byband)
        2390             :       ! Reset band-by-band fluxes
        2391 17584863840 :       if (associated(fluxes%bnd_flux_up)) fluxes%bnd_flux_up(:,:,:) = 0._r8
        2392 17586356112 :       if (associated(fluxes%bnd_flux_dn)) fluxes%bnd_flux_dn(:,:,:) = 0._r8
        2393 17586356112 :       if (associated(fluxes%bnd_flux_net)) fluxes%bnd_flux_net(:,:,:) = 0._r8
        2394 17589340656 :       if (associated(fluxes%bnd_flux_dn_dir)) fluxes%bnd_flux_dn_dir(:,:,:) = 0._r8
        2395             :    end select
        2396             :
        2397     5969088 : end subroutine reset_fluxes
        2398             :
        2399             : !=========================================================================================
        2400             :
        2401     4476816 : subroutine free_optics_sw(optics)
        2402             :
        2403             :    type(ty_optical_props_2str), intent(inout) :: optics
        2404             :
        2405     4476816 :    if (allocated(optics%tau)) deallocate(optics%tau)
        2406     4476816 :    if (allocated(optics%ssa)) deallocate(optics%ssa)
        2407     4476816 :    if (allocated(optics%g)) deallocate(optics%g)
        2408     4476816 :    call optics%finalize()
        2409             :
        2410     4476816 : end subroutine free_optics_sw
        2411             :
        2412             : !=========================================================================================
        2413             :
        2414     2984544 : subroutine free_optics_lw(optics)
        2415             :
        2416             :    type(ty_optical_props_1scl), intent(inout) :: optics
        2417             :
        2418     2984544 :    if (allocated(optics%tau)) deallocate(optics%tau)
        2419     2984544 :    call optics%finalize()
        2420             :
        2421     2984544 : end subroutine free_optics_lw
        2422             :
        2423             : !=========================================================================================
        2424             :
        2425     5969088 : subroutine free_fluxes(fluxes)
        2426             :
        2427             :    class(ty_fluxes_broadband), intent(inout) :: fluxes
        2428             :
        2429     5969088 :    if (associated(fluxes%flux_up)) deallocate(fluxes%flux_up)
        2430     5969088 :    if (associated(fluxes%flux_dn)) deallocate(fluxes%flux_dn)
        2431     5969088 :    if (associated(fluxes%flux_net)) deallocate(fluxes%flux_net)
        2432     5969088 :    if (associated(fluxes%flux_dn_dir)) deallocate(fluxes%flux_dn_dir)
        2433             :
        2434             :    select type (fluxes)
        2435             :    type is (ty_fluxes_byband)
        2436     1492272 :       if (associated(fluxes%bnd_flux_up)) deallocate(fluxes%bnd_flux_up)
        2437     2984544 :       if (associated(fluxes%bnd_flux_dn)) deallocate(fluxes%bnd_flux_dn)
        2438     2984544 :       if (associated(fluxes%bnd_flux_net)) deallocate(fluxes%bnd_flux_net)
        2439     5969088 :       if (associated(fluxes%bnd_flux_dn_dir)) deallocate(fluxes%bnd_flux_dn_dir)
        2440             :    end select
        2441             :
        2442     5969088 : end subroutine free_fluxes
        2443             :
        2444             : !=========================================================================================
        2445             :
        2446      746136 : subroutine modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime)
        2447             :
        2448             :    ! Compute modified cloud fraction, cldfprime.
        2449             :    ! 1. initialize as cld
        2450             :    ! 2. modify for snow if cldfsnow is available. use max(cld, cldfsnow)
        2451             :    ! 3. modify for graupel if cldfgrau is available and graupel_in_rad is true.
        2452             :    !    use max(cldfprime, cldfgrau)
        2453             :
        2454             :    ! Arguments
        2455             :    integer,  intent(in)  :: ncol
        2456             :    real(r8), pointer     :: cld(:,:)       ! cloud fraction
        2457             :    real(r8), pointer     :: cldfsnow(:,:)  ! cloud fraction of just "snow clouds"
        2458             :    real(r8), pointer     :: cldfgrau(:,:)  ! cloud fraction of just "graupel clouds"
        2459             :    real(r8), intent(out) :: cldfprime(:,:) ! modified cloud fraction
        2460             :
        2461             :    ! Local variables
        2462             :    integer :: i, k
        2463             :    !----------------------------------------------------------------------------
        2464             :
        2465      746136 :    if (associated(cldfsnow)) then
        2466    70136784 :       do k = 1, pver
        2467  1159408584 :          do i = 1, ncol
        2468  1158662448 :             cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k))
        2469             :          end do
        2470             :       end do
        2471             :    else
        2472           0 :       cldfprime(:ncol,:) = cld(:ncol,:)
        2473             :    end if
        2474             :
        2475      746136 :    if (associated(cldfgrau) .and. graupel_in_rad) then
        2476           0 :       do k = 1, pver
        2477           0 :          do i = 1, ncol
        2478           0 :             cldfprime(i,k) = max(cldfprime(i,k), cldfgrau(i,k))
        2479             :          end do
        2480             :       end do
        2481             :    end if
        2482             :
        2483      746136 : end subroutine modified_cloud_fraction
        2484             :
        2485             : !=========================================================================================
        2486             :
        2487     9833936 : subroutine stop_on_err(errmsg, sub, info)
        2488             :
        2489             : ! call endrun if RRTMGP function returns non-empty error message.
        2490             :
        2491             :    character(len=*), intent(in) :: errmsg    ! return message from RRTMGP function
        2492             :    character(len=*), intent(in) :: sub       ! name of calling subroutine
        2493             :    character(len=*), intent(in) :: info      ! name of called function
        2494             :
        2495     9833936 :    if (len_trim(errmsg) > 0) then
        2496           0 :       call endrun(trim(sub)//': ERROR: '//trim(info)//': '//trim(errmsg))
        2497             :    end if
        2498             :
        2499     9833936 : end subroutine stop_on_err
        2500             :
        2501             : !=========================================================================================
        2502             :
        2503    17907264 : end module radiation
        2504             :

Generated by: LCOV version 1.14