LCOV - code coverage report
Current view: | top level - physics/camrt - radiation.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 436 | 540 | 80.7 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 13 | 14 | 92.9 %

              Line data    Source code
           1             : module radiation
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : !
           5             : ! CAM interface to the legacy 'camrt' radiation code
           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 physics_types,       only: physics_state, physics_ptend
          13             : use phys_grid,           only: get_ncols_p
          14             : use camsrfexch,          only: cam_out_t, cam_in_t
          15             : use physconst,           only: cpair, cappa
          16             : use time_manager,        only: get_nstep, is_first_restart_step, &
          17             :                                get_curr_calday, get_step_size
          18             : use cam_control_mod,     only: lambm0, obliqr, mvelpp, eccen
          19             :
          20             : use radae,               only: abstot_3d, absnxt_3d, emstot_3d, initialize_radbuffer, ntoplw
          21             :
          22             : use scamMod,             only: scm_crm_mode, single_column,have_cld,cldobs,&
          23             :                                have_clwp,clwpobs,have_tg,tground
          24             :
          25             : use cam_grid_support,    only: cam_grid_write_attr, cam_grid_id,            &
          26             :                                cam_grid_header_info_t, cam_grid_dimensions, &
          27             :                                cam_grid_write_dist_array, cam_grid_read_dist_array
          28             :
          29             : use cam_history,         only: outfld, hist_fld_active
          30             : use cam_history_support, only: fillvalue
          31             :
          32             : use cam_pio_utils,       only: cam_pio_def_dim
          33             :
          34             : use pio,                 only: file_desc_t, var_desc_t, &
          35             :                                pio_double, pio_int, pio_noerr, &
          36             :                                pio_seterrorhandling, pio_bcast_error, &
          37             :                                pio_inq_varid, &
          38             :                                pio_def_var, pio_def_dim, &
          39             :                                pio_put_var, pio_get_var, pio_put_att
          40             :
          41             : use cam_abortutils,      only: endrun
          42             : use error_messages,      only: handle_err
          43             : use perf_mod,            only: t_startf, t_stopf
          44             : use cam_logfile,         only: iulog
          45             :
          46             : implicit none
          47             : private
          48             : save
          49             :
          50             : public :: &
          51             :    radiation_readnl,         &! read namelist variables
          52             :    radiation_register,       &! registers radiation physics buffer fields
          53             :    radiation_nextsw_cday,    &! calendar day of next radiation calculation
          54             :    radiation_do,             &! query which radiation calcs are done this timestep
          55             :    radiation_init,           &! initialization
          56             :    radiation_define_restart, &!
          57             :    radiation_write_restart,  &!
          58             :    radiation_read_restart,   &!
          59             :    radiation_tend,           &! compute heating rates and fluxes
          60             :    rad_out_t                  ! type for diagnostic outputs
          61             :
          62             : real(r8), public, protected :: nextsw_cday       ! future radiation calday for surface models
          63             :
          64             : type rad_out_t
          65             :    real(r8) :: solin(pcols)         ! Solar incident flux
          66             :    real(r8) :: fsntoa(pcols)        ! Net solar flux at TOA
          67             :    real(r8) :: fsutoa(pcols)        ! upwelling solar flux at TOA
          68             :    real(r8) :: fsntoac(pcols)       ! Clear sky net solar flux at TOA
          69             :    real(r8) :: fsnirt(pcols)        ! Near-IR flux absorbed at toa
          70             :    real(r8) :: fsnrtc(pcols)        ! Clear sky near-IR flux absorbed at toa
          71             :    real(r8) :: fsnirtsq(pcols)      ! Near-IR flux absorbed at toa >= 0.7 microns
          72             :    real(r8) :: fsntc(pcols)         ! Clear sky total column abs solar flux
          73             :    real(r8) :: fsnsc(pcols)         ! Clear sky surface abs solar flux
          74             :    real(r8) :: fsdsc(pcols)         ! Clear sky surface downwelling solar flux
          75             :    real(r8) :: flut(pcols)          ! Upward flux at top of model
          76             :    real(r8) :: flutc(pcols)         ! Upward Clear Sky flux at top of model
          77             :    real(r8) :: flntc(pcols)         ! Clear sky lw flux at model top
          78             :    real(r8) :: flnsc(pcols)         ! Clear sky lw flux at srf (up-down)
          79             :    real(r8) :: fldsc(pcols)         ! Clear sky lw flux at srf (down)
          80             :    real(r8) :: flwds(pcols)         ! Down longwave flux at surface
          81             :    real(r8) :: fsnr(pcols)
          82             :    real(r8) :: flnr(pcols)
          83             :    real(r8) :: fsds(pcols)          ! Surface solar down flux
          84             :    real(r8) :: fln200(pcols)        ! net longwave flux interpolated to 200 mb
          85             :    real(r8) :: fln200c(pcols)       ! net clearsky longwave flux interpolated to 200 mb
          86             :    real(r8) :: fsn200(pcols)        ! fns interpolated to 200 mb
          87             :    real(r8) :: fsn200c(pcols)       ! fcns interpolated to 200 mb
          88             :    real(r8) :: sols(pcols)          ! Solar downward visible direct  to surface
          89             :    real(r8) :: soll(pcols)          ! Solar downward near infrared direct  to surface
          90             :    real(r8) :: solsd(pcols)         ! Solar downward visible diffuse to surface
          91             :    real(r8) :: solld(pcols)         ! Solar downward near infrared diffuse to surface
          92             :    real(r8) :: qrsc(pcols,pver)     ! clearsky shortwave radiative heating rate
          93             :    real(r8) :: qrlc(pcols,pver)     ! clearsky longwave  radiative heating rate
          94             :    real(r8) :: fsdtoa(pcols)        ! Solar input = Flux Solar Downward Top of Atmosphere
          95             :    real(r8) :: swcf(pcols)          ! shortwave cloud forcing
          96             :    real(r8) :: lwcf(pcols)          ! longwave cloud forcing
          97             :
          98             :    real(r8) :: tot_cld_vistau(pcols,pver)
          99             :    real(r8) :: tot_icld_vistau(pcols,pver)
         100             :    real(r8) :: liq_icld_vistau(pcols,pver)  ! in-cld liq cloud optical depth (only during day, night = fillvalue)
         101             :    real(r8) :: ice_icld_vistau(pcols,pver)  ! in-cld ice cloud optical depth (only during day, night = fillvalue)
         102             : end type rad_out_t
         103             :
         104             : ! Namelist variables
         105             :
         106             : character(len=cl) :: absems_data
         107             : integer :: iradsw = -1     ! freq. of shortwave radiation calc in time steps (positive)
         108             :                            ! or hours (negative).
         109             : integer :: iradlw = -1     ! frequency of longwave rad. calc. in time steps (positive)
         110             :                            ! or hours (negative).
         111             : integer :: iradae = -12    ! frequency of absorp/emis calc in time steps (positive)
         112             :                            ! or hours (negative).
         113             : integer :: irad_always = 0 ! Specifies length of time in timesteps (positive)
         114             :                            ! or hours (negative) SW/LW radiation will be
         115             :                            ! run continuously from the start of an
         116             :                            ! initial or restart run
         117             : logical :: use_rad_dt_cosz = .false. ! if true use zenith angle averaged over
         118             :                                      ! interval between radiation calculations
         119             :
         120             : ! Physics buffer indices
         121             : integer :: qrs_idx      = 0
         122             : integer :: qrl_idx      = 0
         123             : integer :: fsds_idx     = 0
         124             : integer :: fsns_idx     = 0
         125             : integer :: fsnt_idx     = 0
         126             : integer :: flns_idx     = 0
         127             : integer :: flnt_idx     = 0
         128             : integer :: cld_idx      = 0
         129             : integer :: rel_idx      = 0
         130             : integer :: rei_idx      = 0
         131             : integer :: cicewp_idx   = -1
         132             : integer :: cliqwp_idx   = -1
         133             : integer :: cldemis_idx  = -1
         134             : integer :: cldtau_idx   = -1
         135             : integer :: nmxrgn_idx   = -1
         136             : integer :: pmxrgn_idx   = -1
         137             :
         138             : ! averaging time interval for zenith angle
         139             : real(r8) :: dt_avg = 0._r8
         140             :
         141             : real(r8), parameter :: cgs2mks = 1.e-3_r8
         142             :
         143             : ! PIO descriptors (for restarts)
         144             :
         145             : type(var_desc_t), allocatable :: abstot_desc(:)
         146             : type(var_desc_t) :: emstot_desc, absnxt_desc(4)
         147             : type(var_desc_t) :: nextsw_cday_desc
         148             :
         149             : logical  :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the zenith calculation
         150             : real(r8) :: rad_uniform_angle = -99._r8
         151             :
         152             : !===============================================================================
         153             : contains
         154             : !===============================================================================
         155             :
         156        1536 : subroutine radiation_readnl(nlfile)
         157             :
         158             :    ! Read radiation_nl namelist group.
         159             :
         160             :    use namelist_utils,  only: find_group_name
         161             :    use units,           only: getunit, freeunit
         162             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, &
         163             :                               mpi_character, mpi_real8
         164             :
         165             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         166             :
         167             :    ! Local variables
         168             :    integer :: unitn, ierr
         169             :    integer :: dtime      ! timestep size
         170             :    character(len=*), parameter :: sub = 'radiation_readnl'
         171             :
         172             :    namelist /radiation_nl/ absems_data, iradsw, iradlw, iradae, irad_always, &
         173             :                            use_rad_dt_cosz, use_rad_uniform_angle, rad_uniform_angle
         174             :    !-----------------------------------------------------------------------------
         175             :
         176        1536 :    if (masterproc) then
         177           2 :       unitn = getunit()
         178           2 :       open( unitn, file=trim(nlfile), status='old' )
         179           2 :       call find_group_name(unitn, 'radiation_nl', status=ierr)
         180           2 :       if (ierr == 0) then
         181           2 :          read(unitn, radiation_nl, iostat=ierr)
         182           2 :          if (ierr /= 0) then
         183           0 :             call endrun(sub // ':: ERROR reading namelist')
         184             :          end if
         185             :       end if
         186           2 :       close(unitn)
         187           2 :       call freeunit(unitn)
         188             :    end if
         189             :
         190             :    ! Broadcast namelist variables
         191        1536 :    call mpi_bcast(absems_data, len(absems_data), mpi_character, mstrid, mpicom, ierr)
         192        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: absems_data")
         193        1536 :    call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr)
         194        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw")
         195        1536 :    call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr)
         196        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw")
         197        1536 :    call mpi_bcast(iradae, 1, mpi_integer, mstrid, mpicom, ierr)
         198        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradae")
         199        1536 :    call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr)
         200        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always")
         201        1536 :    call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr)
         202        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz")
         203        1536 :    call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr)
         204        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle")
         205        1536 :    call mpi_bcast(rad_uniform_angle, 1, mpi_real8,  mstrid, mpicom, ierr)
         206        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle")
         207             :
         208        1536 :    if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then
         209           0 :       call endrun(sub // ' ERROR - use_rad_uniform_angle is set to .true, but rad_uniform_angle is not set ')
         210             :    end if
         211             :
         212             :    ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
         213        1536 :    dtime  = get_step_size()
         214        1536 :    if (iradsw      < 0) iradsw      = nint((-iradsw     *3600._r8)/dtime)
         215        1536 :    if (iradlw      < 0) iradlw      = nint((-iradlw     *3600._r8)/dtime)
         216        1536 :    if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime)
         217             :
         218             :    ! Convert iradae from hours to timesteps if necessary and check that
         219             :    ! iradae must be an even multiple of iradlw
         220        1536 :    if (iradae < 0) iradae = nint((-iradae*3600._r8)/dtime)
         221        1536 :    if (mod(iradae,iradlw)/=0) then
         222           0 :       write(iulog,*) sub//': iradae must be an even multiple of iradlw.'
         223           0 :       write(iulog,*)'     iradae = ',iradae,', iradlw = ',iradlw
         224           0 :       call endrun(sub//': iradae must be an even multiple of iradlw.')
         225             :    end if
         226             :
         227             :    !-----------------------------------------------------------------------
         228             :    ! Print runtime options to log.
         229             :    !-----------------------------------------------------------------------
         230             :
         231        1536 :    if (masterproc) then
         232           2 :       write(iulog,*) 'CAMRT radiation scheme parameters:'
         233           2 :       write(iulog,10) iradsw, iradlw, iradae, irad_always, use_rad_dt_cosz
         234           2 :       write(iulog,*) '  Abs/Emis dataset: ', trim(absems_data)
         235             :    end if
         236             :
         237             : 10 format('  Frequency (timesteps) of Shortwave Radiation calc:    ',i5/, &
         238             :           '  Frequency (timesteps) of Longwave Radiation calc:     ',i5/, &
         239             :           '  Frequency (timesteps) of Absorptivity/Emissivity calc:',i5/, &
         240             :           '  SW/LW calc done every timestep for first N steps. N=  ',i5/, &
         241             :           '  Use average zenith angle:                             ',l5)
         242             :
         243             :
         244        1536 : end subroutine radiation_readnl
         245             :
         246             : !================================================================================================
         247             :
         248        1536 : subroutine radiation_register
         249             :
         250             :    ! Register radiation fields in the physics buffer
         251             :
         252             :    use physics_buffer, only: pbuf_add_field, dtype_r8
         253             :    use radiation_data, only: rad_data_register
         254             :
         255        1536 :    call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate
         256        1536 :    call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave  radiative heating rate
         257             :
         258        1536 :    call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux
         259             :
         260        1536 :    call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux
         261        1536 :    call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux
         262        1536 :    call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux
         263        1536 :    call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux
         264             :
         265        1536 :    call rad_data_register()
         266             :
         267        1536 : end subroutine radiation_register
         268             :
         269             : !================================================================================================
         270             :
         271     9731472 : function radiation_do(op, timestep)
         272             :
         273             :    ! Returns true if the specified operation is done this timestep.
         274             :
         275             :    character(len=*), intent(in) :: op             ! name of operation
         276             :    integer, intent(in), optional:: timestep
         277             :    logical                      :: radiation_do   ! return value
         278             :
         279             :    ! Local variables
         280             :    integer :: nstep             ! current timestep number
         281             :    !-----------------------------------------------------------------------
         282             :
         283     9731472 :    if (present(timestep)) then
         284     2250792 :       nstep = timestep
         285             :    else
         286     7480680 :       nstep = get_nstep()
         287             :    end if
         288             :
         289     5241528 :    select case (op)
         290             :
         291             :    case ('sw') ! do a shortwave heating calc this timestep?
         292             :       radiation_do = nstep == 0  .or.  iradsw == 1                     &
         293             :                     .or. (mod(nstep-1,iradsw) == 0  .and.  nstep /= 1) &
         294     5241528 :                     .or. nstep <= irad_always
         295             :
         296             :    case ('lw') ! do a longwave heating calc this timestep?
         297             :       radiation_do = nstep == 0  .or.  iradlw == 1                     &
         298             :                     .or. (mod(nstep-1,iradlw) == 0  .and.  nstep /= 1) &
         299     2990736 :                     .or. nstep <= irad_always
         300             :
         301             :    case ('absems') ! do an absorptivity/emissivity calculation this timestep?
         302             :       radiation_do = nstep == 0  .or.  iradae == 1                     &
         303     1495368 :                     .or. (mod(nstep-1,iradae) == 0  .and.  nstep /= 1)
         304             :
         305             :    case ('aeres') ! write absorptivity/emissivity to restart file this timestep?
         306        3840 :       radiation_do = mod(nstep,iradae) /= 0
         307             :
         308             :    case default
         309     9731472 :       call endrun('radiation_do: unknown operation:'//op)
         310             :
         311             :    end select
         312        1536 : end function radiation_do
         313             :
         314             : !================================================================================================
         315             :
         316     1495368 : real(r8) function radiation_nextsw_cday()
         317             :
         318             :    ! Returns calendar day of next sw radiation calculation
         319             :
         320             :    ! Local variables
         321             :    integer :: nstep      ! timestep counter
         322             :    logical :: dosw       ! true => do shosrtwave calc
         323             :    integer :: offset     ! offset for calendar day calculation
         324             :    integer :: dtime      ! integer timestep size
         325             :    real(r8):: calday     ! calendar day of
         326             :    real(r8):: caldayp1   ! calendar day of next time-step
         327             :    !-----------------------------------------------------------------------
         328             :
         329     1495368 :    radiation_nextsw_cday = -1._r8
         330     1495368 :    dosw   = .false.
         331     1495368 :    nstep  = get_nstep()
         332     1495368 :    dtime  = get_step_size()
         333     1495368 :    offset = 0
         334             :    do while (.not. dosw)
         335     2250792 :       nstep = nstep + 1
         336     2250792 :       offset = offset + dtime
         337     2250792 :       if (radiation_do('sw', nstep)) then
         338     1495368 :          radiation_nextsw_cday = get_curr_calday(offset=offset)
         339             :          dosw = .true.
         340             :       end if
         341             :    end do
         342     1495368 :    if(radiation_nextsw_cday == -1._r8) then
         343           0 :       call endrun('error in radiation_nextsw_cday')
         344             :    end if
         345             :
         346             :    ! determine if next radiation time-step not equal to next time-step
         347     1495368 :    if (get_nstep() >= 1) then
         348     1492272 :       caldayp1 = get_curr_calday(offset=int(dtime))
         349     1492272 :       if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8
         350             :    end if
         351             :
         352     1495368 : end function radiation_nextsw_cday
         353             :
         354             : !================================================================================================
         355             :
         356        1536 : subroutine radiation_init(pbuf2d)
         357             :
         358             :    ! Initialize the radiation parameterization, add fields to the history buffer
         359             :
         360             :    use cam_history,    only: addfld, add_default, horiz_only
         361             :    use physconst,      only: gravit, cpair, epsilo, stebol, &
         362             :                              pstd, mwdry, mwco2, mwo3
         363             :
         364             :    use physics_buffer, only: physics_buffer_desc, pbuf_get_index
         365             :    use radsw,          only: radsw_init
         366             :    use radlw,          only: radlw_init
         367             :    use radae,          only: radae_init
         368             :    use radconstants,   only: radconstants_init
         369             :    use rad_solar_var,  only: rad_solar_var_init
         370             :    use radiation_data, only: rad_data_init
         371             :    use phys_control,   only: phys_getopts
         372             :    use time_manager,   only: get_step_size
         373             :
         374             :    ! args
         375             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         376             :
         377             :
         378             :    ! Local variables
         379             :    integer :: nstep                       ! current timestep number
         380             :    logical :: history_amwg                ! output the variables used by the AMWG diag package
         381             :    logical :: history_vdiag               ! output the variables used by the AMWG variability diag package
         382             :    logical :: history_budget              ! output tendencies and state variables for CAM4
         383             :                                           ! temperature, water vapor, cloud ice and cloud
         384             :                                           ! liquid budgets.
         385             :    integer :: history_budget_histfile_num ! output history file number for budget fields
         386             :
         387             :    integer :: dtime
         388             :    !-----------------------------------------------------------------------
         389             :
         390        1536 :    call radconstants_init()
         391        1536 :    call rad_solar_var_init()
         392             :
         393        1536 :    call radsw_init(gravit)
         394        1536 :    call radlw_init(gravit, stebol)
         395             :    call radae_init( &
         396             :       gravit, epsilo, stebol, pstd, mwdry, &
         397        1536 :       mwco2, mwo3, absems_data)
         398             :
         399        1536 :    call rad_data_init(pbuf2d)
         400             :
         401             :    ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation
         402        1536 :    if (use_rad_dt_cosz)  then
         403           0 :       dtime  = get_step_size()
         404           0 :       dt_avg = real(iradsw*dtime, r8)
         405             :    end if
         406             :
         407             :    ! Surface components to get radiation computed today
         408        1536 :    if (.not. is_first_restart_step()) then
         409         768 :       nextsw_cday = get_curr_calday()
         410             :    end if
         411             :
         412             :    ! Get physics buffer indices
         413        1536 :    cld_idx    = pbuf_get_index('CLD')
         414        1536 :    rel_idx    = pbuf_get_index('REL')
         415        1536 :    rei_idx    = pbuf_get_index('REI')
         416             :
         417             :    ! "irad_always" is number of time steps to execute radiation continuously from start of
         418             :    ! initial OR restart run
         419        1536 :    nstep = get_nstep()
         420        1536 :    if ( irad_always > 0) then
         421           0 :       nstep       = get_nstep()
         422           0 :       irad_always = irad_always + nstep
         423             :    end if
         424             :
         425             :    ! Shortwave radiation
         426        1536 :    call addfld ('SOLIN',           horiz_only,   'A','W/m2','Solar insolation',                            sampling_seq='rad_lwsw')
         427             :    call addfld ('SOLL',            horiz_only,   'A','W/m2','Solar downward near infrared direct  to surface',              &
         428        1536 :                                                                                                            sampling_seq='rad_lwsw')
         429        1536 :    call addfld ('SOLS',            horiz_only,   'A','W/m2','Solar downward visible direct  to surface',   sampling_seq='rad_lwsw')
         430             :    call addfld ('SOLLD',           horiz_only,   'A','W/m2','Solar downward near infrared diffuse to surface',              &
         431        1536 :                                                                                                            sampling_seq='rad_lwsw')
         432        1536 :    call addfld ('SOLSD',           horiz_only,   'A','W/m2','Solar downward visible diffuse to surface',   sampling_seq='rad_lwsw')
         433        3072 :    call addfld ('QRS',             (/ 'lev' /),  'A','K/s', 'Solar heating rate',                          sampling_seq='rad_lwsw')
         434        3072 :    call addfld ('QRSC',            (/ 'lev' /),  'A','K/s', 'Clearsky solar heating rate',                 sampling_seq='rad_lwsw')
         435        1536 :    call addfld ('FSNS',            horiz_only,   'A','W/m2','Net solar flux at surface',                   sampling_seq='rad_lwsw')
         436        1536 :    call addfld ('FSNT',            horiz_only,   'A','W/m2','Net solar flux at top of model',              sampling_seq='rad_lwsw')
         437        1536 :    call addfld ('FSNTOA',          horiz_only,   'A','W/m2','Net solar flux at top of atmosphere',         sampling_seq='rad_lwsw')
         438        1536 :    call addfld ('FSUTOA',          horiz_only,   'A','W/m2','Upwelling solar flux at top of atmosphere',   sampling_seq='rad_lwsw')
         439             :    call addfld ('FSNTOAC',         horiz_only,   'A','W/m2','Clearsky net solar flux at top of atmosphere',                 &
         440        1536 :                                                                                                            sampling_seq='rad_lwsw')
         441        1536 :    call addfld ('FSDTOA',          horiz_only,   'A','W/m2','Downwelling solar flux at top of atmosphere', sampling_seq='rad_lwsw')
         442        1536 :    call addfld ('FSN200',          horiz_only,   'A','W/m2','Net shortwave flux at 200 mb',                sampling_seq='rad_lwsw')
         443        1536 :    call addfld ('FSN200C',         horiz_only,   'A','W/m2','Clearsky net shortwave flux at 200 mb',       sampling_seq='rad_lwsw')
         444        1536 :    call addfld ('FSNTC',           horiz_only,   'A','W/m2','Clearsky net solar flux at top of model',     sampling_seq='rad_lwsw')
         445        1536 :    call addfld ('FSNSC',           horiz_only,   'A','W/m2','Clearsky net solar flux at surface',          sampling_seq='rad_lwsw')
         446             :    call addfld ('FSDSC',           horiz_only,   'A','W/m2','Clearsky downwelling solar flux at surface',                   &
         447        1536 :                                                                                                            sampling_seq='rad_lwsw')
         448        1536 :    call addfld ('FSDS',            horiz_only,   'A','W/m2','Downwelling solar flux at surface',           sampling_seq='rad_lwsw')
         449        3072 :    call addfld ('FUS',             (/ 'ilev' /), 'I','W/m2','Shortwave upward flux')
         450        3072 :    call addfld ('FDS',             (/ 'ilev' /), 'I','W/m2','Shortwave downward flux')
         451        3072 :    call addfld ('FUSC',            (/ 'ilev' /), 'I','W/m2','Shortwave clear-sky upward flux')
         452        3072 :    call addfld ('FDSC',            (/ 'ilev' /), 'I','W/m2','Shortwave clear-sky downward flux')
         453             :    call addfld ('FSNIRTOA',        horiz_only,   'A','W/m2','Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere',  &
         454        1536 :                                                                                                            sampling_seq='rad_lwsw')
         455             :    call addfld ('FSNRTOAC',        horiz_only,   'A','W/m2',                                                                &
         456        1536 :          'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere',                           sampling_seq='rad_lwsw')
         457             :    call addfld ('FSNRTOAS',        horiz_only,   'A','W/m2','Net near-infrared flux (>= 0.7 microns) at top of atmosphere', &
         458        1536 :                                                                                                            sampling_seq='rad_lwsw')
         459        1536 :    call addfld ('FSNR',            horiz_only,   'A','W/m2','Net solar flux at tropopause',                sampling_seq='rad_lwsw')
         460        1536 :    call addfld ('SWCF',            horiz_only,   'A','W/m2','Shortwave cloud forcing',                     sampling_seq='rad_lwsw')
         461             :
         462             :    call addfld ('TOT_CLD_VISTAU',  (/ 'lev' /),  'A','1',   'Total gbx cloud visible sw optical depth',                     &
         463        3072 :                                                                                          sampling_seq='rad_lwsw',flag_xyfill=.true.)
         464             :    call addfld ('TOT_ICLD_VISTAU', (/ 'lev' /),  'A','1',   'Total in-cloud visible sw optical depth',                      &
         465        3072 :                                                                                          sampling_seq='rad_lwsw',flag_xyfill=.true.)
         466             :    call addfld ('LIQ_ICLD_VISTAU', (/ 'lev' /),  'A','1',   'Liquid in-cloud visible sw optical depth',                     &
         467        3072 :                                                                                          sampling_seq='rad_lwsw',flag_xyfill=.true.)
         468             :    call addfld ('ICE_ICLD_VISTAU', (/ 'lev' /),  'A','1',   'Ice in-cloud visible sw optical depth',                        &
         469        3072 :                                                                                          sampling_seq='rad_lwsw',flag_xyfill=.true.)
         470             :
         471             :    ! Longwave radiation
         472        3072 :    call addfld ('QRL',             (/ 'lev' /),  'A','K/s', 'Longwave heating rate',                       sampling_seq='rad_lwsw')
         473        3072 :    call addfld ('QRLC',            (/ 'lev' /),  'A','K/s', 'Clearsky longwave heating rate',              sampling_seq='rad_lwsw')
         474        1536 :    call addfld ('FLNS',            horiz_only,   'A','W/m2','Net longwave flux at surface',                sampling_seq='rad_lwsw')
         475        1536 :    call addfld ('FLDS',            horiz_only,   'A','W/m2','Downwelling longwave flux at surface',        sampling_seq='rad_lwsw')
         476        1536 :    call addfld ('FLNT',            horiz_only,   'A','W/m2','Net longwave flux at top of model',           sampling_seq='rad_lwsw')
         477        1536 :    call addfld ('FLUT',            horiz_only,   'A','W/m2','Upwelling longwave flux at top of model',     sampling_seq='rad_lwsw')
         478             :    call addfld ('FLUTC',           horiz_only,   'A','W/m2','Clearsky upwelling longwave flux at top of model',             &
         479        1536 :                                                                                                            sampling_seq='rad_lwsw')
         480        1536 :    call addfld ('FLNTC',           horiz_only,   'A','W/m2','Clearsky net longwave flux at top of model',  sampling_seq='rad_lwsw')
         481        1536 :    call addfld ('FLN200',          horiz_only,   'A','W/m2','Net longwave flux at 200 mb',                 sampling_seq='rad_lwsw')
         482        1536 :    call addfld ('FLN200C',         horiz_only,   'A','W/m2','Clearsky net longwave flux at 200 mb',        sampling_seq='rad_lwsw')
         483        1536 :    call addfld ('FLNR',            horiz_only,   'A','W/m2','Net longwave flux at tropopause',             sampling_seq='rad_lwsw')
         484        1536 :    call addfld ('FLNSC',           horiz_only,   'A','W/m2','Clearsky net longwave flux at surface',       sampling_seq='rad_lwsw')
         485             :    call addfld ('FLDSC',           horiz_only,   'A','W/m2','Clearsky downwelling longwave flux at surface',                &
         486        1536 :                                                                                                            sampling_seq='rad_lwsw')
         487        1536 :    call addfld ('LWCF',            horiz_only,   'A','W/m2','Longwave cloud forcing',                      sampling_seq='rad_lwsw')
         488        3072 :    call addfld ('FUL',             (/ 'ilev' /), 'I','W/m2','Longwave upward flux')
         489        3072 :    call addfld ('FDL',             (/ 'ilev' /), 'I','W/m2','Longwave downward flux')
         490        3072 :    call addfld ('FULC',            (/ 'ilev' /), 'I','W/m2','Longwave clear-sky upward flux')
         491        3072 :    call addfld ('FDLC',            (/ 'ilev' /), 'I','W/m2','Longwave clear-sky downward flux')
         492             :
         493             :    ! Heating rate needed for d(theta)/dt computation
         494        3072 :    call addfld ('HR',              (/ 'lev' /),  'A','K/s', 'Heating rate needed for d(theta)/dt computation')
         495             :
         496             :    ! determine default variables
         497             :    call phys_getopts(history_amwg_out   = history_amwg,   &
         498             :                       history_vdiag_out  = history_vdiag,  &
         499             :                       history_budget_out = history_budget, &
         500        1536 :                       history_budget_histfile_num_out = history_budget_histfile_num)
         501             :
         502        1536 :    if (history_amwg) then
         503             :       ! Shortwave variables
         504        1536 :       call add_default ('SOLIN   ', 1, ' ')
         505        1536 :       call add_default ('QRS     ', 1, ' ')
         506        1536 :       call add_default ('FSNS    ', 1, ' ')
         507        1536 :       call add_default ('FSNT    ', 1, ' ')
         508        1536 :       call add_default ('FSDTOA  ', 1, ' ')
         509        1536 :       call add_default ('FSNTOA  ', 1, ' ')
         510        1536 :       call add_default ('FSUTOA  ', 1, ' ')
         511        1536 :       call add_default ('FSNTOAC ', 1, ' ')
         512        1536 :       call add_default ('FSNTC   ', 1, ' ')
         513        1536 :       call add_default ('FSNSC   ', 1, ' ')
         514        1536 :       call add_default ('FSDSC   ', 1, ' ')
         515        1536 :       call add_default ('FSDS    ', 1, ' ')
         516        1536 :       call add_default ('SWCF    ', 1, ' ')
         517             :       ! Longwave variables
         518        1536 :       call add_default ('QRL     ', 1, ' ')
         519        1536 :       call add_default ('FLNS    ', 1, ' ')
         520        1536 :       call add_default ('FLDS    ', 1, ' ')
         521        1536 :       call add_default ('FLNT    ', 1, ' ')
         522        1536 :       call add_default ('FLUT    ', 1, ' ')
         523        1536 :       call add_default ('FLUTC   ', 1, ' ')
         524        1536 :       call add_default ('FLNTC   ', 1, ' ')
         525        1536 :       call add_default ('FLNSC   ', 1, ' ')
         526        1536 :       call add_default ('FLDSC   ', 1, ' ')
         527        1536 :       call add_default ('LWCF    ', 1, ' ')
         528             :    endif
         529        1536 :    if (single_column.and.scm_crm_mode) then
         530             :       ! Shortwave variables
         531           0 :       call add_default ('FUS     ', 1, ' ')
         532           0 :       call add_default ('FUSC    ', 1, ' ')
         533           0 :       call add_default ('FDS     ', 1, ' ')
         534           0 :       call add_default ('FDSC    ', 1, ' ')
         535             :       ! Longwave variables
         536           0 :       call add_default ('FUL     ', 1, ' ')
         537           0 :       call add_default ('FULC    ', 1, ' ')
         538           0 :       call add_default ('FDL     ', 1, ' ')
         539           0 :       call add_default ('FDLC    ', 1, ' ')
         540             :    endif
         541             :
         542        1536 :    if ( history_budget .and. history_budget_histfile_num > 1 ) then
         543           0 :       call add_default ('QRL     ', history_budget_histfile_num, ' ')
         544           0 :       call add_default ('QRS     ', history_budget_histfile_num, ' ')
         545             :    end if
         546             :
         547        1536 :    if (history_vdiag) then
         548           0 :       call add_default('FLUT',2,' ')
         549           0 :       call add_default('FLUT',3,' ')
         550             :    end if
         551             :
         552        1536 :    cicewp_idx = pbuf_get_index('CICEWP')
         553        1536 :    cliqwp_idx = pbuf_get_index('CLIQWP')
         554        1536 :    cldemis_idx= pbuf_get_index('CLDEMIS')
         555        1536 :    cldtau_idx = pbuf_get_index('CLDTAU')
         556        1536 :    nmxrgn_idx = pbuf_get_index('NMXRGN')
         557        1536 :    pmxrgn_idx = pbuf_get_index('PMXRGN')
         558             :
         559        1536 : end subroutine radiation_init
         560             :
         561             : !===============================================================================
         562             :
         563        1536 : subroutine radiation_define_restart(file)
         564             :
         565             :    ! define variables to be written to restart file
         566             :
         567             :    ! arguments
         568             :    type(file_desc_t), intent(inout) :: file
         569             :
         570             :    ! local variables
         571             :    integer :: i, ierr
         572             :    integer :: grid_id
         573             :    integer :: hdimcnt
         574             :    integer :: pver_id, pverp_id
         575             :    integer :: vsize
         576             :    integer :: dimids(4)
         577             :
         578        1536 :    type(cam_grid_header_info_t) :: info
         579             :
         580             :    character(len=16) :: pname
         581             :    !----------------------------------------------------------------------------
         582             :
         583        1536 :    call pio_seterrorhandling(File, PIO_BCAST_ERROR)
         584             :
         585        1536 :    ierr = pio_def_var(File, 'nextsw_cday', pio_double, nextsw_cday_desc)
         586        1536 :    ierr = pio_put_att(File, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models')
         587             :
         588        1536 :    if (radiation_do('aeres')) then
         589             :
         590           0 :       grid_id = cam_grid_id('physgrid')
         591           0 :       call cam_grid_write_attr(File, grid_id, info)
         592           0 :       hdimcnt = info%num_hdims()
         593           0 :       do i = 1, hdimcnt
         594           0 :          dimids(i) = info%get_hdimid(i)
         595             :       end do
         596             :
         597           0 :       call cam_pio_def_dim(File, 'lev',  pver,  pver_id,  existOK=.true.)
         598           0 :       call cam_pio_def_dim(File, 'ilev', pverp, pverp_id, existOK=.true.)
         599             :
         600           0 :       vsize = pverp - ntoplw + 1
         601           0 :       if (vsize /= pverp) then
         602           0 :          ierr = pio_def_dim(File, 'lwcols', vsize, dimids(hdimcnt+1))
         603             :       else
         604           0 :          dimids(hdimcnt+1) = pverp_id
         605             :       end if
         606             :
         607             :       ! split into vsize variables to avoid excessive memory usage in IO
         608             :
         609           0 :       allocate(abstot_desc(ntoplw:pverp))
         610             :
         611           0 :       do i = ntoplw, pverp
         612           0 :          write(pname,'(a,i3.3)') 'NAL_absorp', i
         613           0 :          ierr = pio_def_var(File, trim(pname), pio_double, dimids(1:hdimcnt+1), abstot_desc(i))
         614             :       end do
         615             :
         616           0 :       dimids(hdimcnt+1) = pverp_id
         617           0 :       ierr = pio_def_var(File, 'Emissivity', pio_double, dimids(1:hdimcnt+1), emstot_desc)
         618             :
         619           0 :       dimids(hdimcnt+1) = pver_id
         620           0 :       do i=1,4
         621           0 :          write(pname,'(a,i3.3)') 'NN_absorp',i
         622           0 :          ierr = pio_def_var(File, pname, pio_double, dimids(1:hdimcnt+1), absnxt_desc(i))
         623             :       end do
         624             :
         625             :    end if
         626             :
         627        3072 : end subroutine radiation_define_restart
         628             :
         629             : !===============================================================================
         630             :
         631        1536 : subroutine radiation_write_restart(file)
         632             :
         633             :    ! write variables to restart file
         634             :
         635             :    ! arguments
         636             :    type(file_desc_t), intent(inout) :: file
         637             :
         638             :    ! local variables
         639             :    integer :: i, ierr
         640             :    integer :: physgrid
         641             :    integer :: dims(3), gdims(3)
         642             :    integer :: nhdims
         643             :    integer :: ncol
         644             :    !----------------------------------------------------------------------------
         645             :
         646        3072 :    ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /))
         647             :
         648        1536 :    if ( radiation_do('aeres')  ) then
         649             :
         650           0 :       physgrid = cam_grid_id('physgrid')
         651           0 :       call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)
         652             :
         653           0 :       do i = begchunk, endchunk
         654           0 :          ncol = get_ncols_p(i)
         655           0 :          if (ncol < pcols) then
         656           0 :             abstot_3d(ncol+1:pcols,:,:,i) = fillvalue
         657           0 :             absnxt_3d(ncol+1:pcols,:,:,i) = fillvalue
         658           0 :             emstot_3d(ncol+1:pcols,:,i)   = fillvalue
         659             :          end if
         660             :       end do
         661             :
         662             :       ! abstot_3d is written as a series of 3D variables
         663             :
         664           0 :       dims(1) = size(abstot_3d, 1) ! Should be pcols
         665           0 :       dims(2) = size(abstot_3d, 2) ! Should be (pverp-ntoplw+1)
         666           0 :       dims(3) = size(abstot_3d, 4) ! Should be endchunk - begchunk + 1
         667           0 :       gdims(nhdims+1) = dims(2)
         668           0 :       do i = ntoplw, pverp
         669             :          call cam_grid_write_dist_array(File, physgrid, dims(1:3),             &
         670           0 :               gdims(1:nhdims+1), abstot_3d(:,:,i,:), abstot_desc(i))
         671             :       end do
         672             :
         673           0 :       dims(1) = size(emstot_3d, 1) ! Should be pcols
         674           0 :       dims(2) = size(emstot_3d, 2) ! Should be pverp
         675           0 :       dims(3) = size(emstot_3d, 3) ! Should be endchunk - begchunk + 1
         676           0 :       gdims(nhdims+1) = dims(2)
         677             :       call cam_grid_write_dist_array(File, physgrid, dims(1:3),               &
         678           0 :            gdims(1:nhdims+1), emstot_3d, emstot_desc)
         679             :
         680           0 :       dims(1) = size(absnxt_3d, 1) ! Should be pcols
         681           0 :       dims(2) = size(absnxt_3d, 2) ! Should be pver
         682           0 :       dims(3) = size(absnxt_3d, 4) ! Should be endchunk - begchunk + 1
         683           0 :       gdims(nhdims+1) = dims(2)
         684           0 :       do i = 1, 4
         685             :          call cam_grid_write_dist_array(File, physgrid, dims(1:3),             &
         686           0 :               gdims(1:nhdims+1), absnxt_3d(:,:,i,:), absnxt_desc(i))
         687             :       end do
         688             :
         689             :       ! module data was allocated in radiation_define_restart
         690           0 :       deallocate(abstot_desc)
         691             :    end if
         692             :
         693        1536 : end subroutine radiation_write_restart
         694             :
         695             : !===============================================================================
         696             :
         697         768 : subroutine radiation_read_restart(file)
         698             :
         699             :    ! read variables from restart file
         700             :
         701             :    ! arguments
         702             :    type(file_desc_t), intent(inout) :: file
         703             :
         704             :    ! local variables
         705             :
         706             :    integer :: err_handling
         707             :    integer :: ierr
         708             :    integer :: physgrid
         709             :    integer :: dims(3), gdims(3), nhdims
         710             :    integer :: vsize
         711             :    integer :: i
         712             :    real(r8) :: temp_var
         713             :
         714             :    type(var_desc_t) :: vardesc
         715             :    character(len=16) :: pname
         716             :    !----------------------------------------------------------------------------
         717             :
         718             :    ! Put this call here for now.  It should move to an init method when the
         719             :    ! initialization and restart sequencing is unified.
         720         768 :    call initialize_radbuffer()
         721             :
         722         768 :    if ( radiation_do('aeres')  ) then
         723             :
         724         768 :       call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
         725         768 :       ierr = pio_inq_varid(File, 'Emissivity', vardesc)
         726         768 :       call pio_seterrorhandling(File, err_handling)
         727         768 :       if (ierr /= PIO_NOERR) then
         728         768 :          if (masterproc) write(iulog,*) 'Warning: Emissivity variable not found on restart file.'
         729         768 :          return
         730             :       end if
         731             :
         732           0 :       physgrid = cam_grid_id('physgrid')
         733           0 :       call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)
         734             :
         735           0 :       dims(1) = pcols
         736           0 :       dims(2) = pverp
         737           0 :       dims(3) = endchunk - begchunk + 1
         738           0 :       gdims(nhdims+1) = dims(2)
         739             :
         740             :       call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
         741           0 :            gdims(1:nhdims+1), emstot_3d, vardesc)
         742             :
         743           0 :       vsize = pverp - ntoplw + 1
         744           0 :       dims(2) = vsize
         745           0 :       gdims(nhdims+1) = dims(2)
         746             :
         747           0 :       do i = ntoplw, pverp
         748           0 :          write(pname,'(a,i3.3)') 'NAL_absorp', i
         749           0 :          ierr = pio_inq_varid(File, trim(pname), vardesc)
         750             :          call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
         751           0 :               gdims(1:nhdims+1), abstot_3d(:,:,i,:), vardesc)
         752             :       end do
         753             :
         754           0 :       dims(2) = pver
         755           0 :       gdims(nhdims+1) = dims(2)
         756           0 :       do i = 1, 4
         757           0 :          write(pname,'(a,i3.3)') 'NN_absorp', i
         758           0 :          ierr = pio_inq_varid(File, trim(pname), vardesc)
         759             :          call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
         760           0 :               gdims(1:nhdims+1), absnxt_3d(:,:,i,:), vardesc)
         761             :       end do
         762             :    end if
         763             :
         764           0 :    ierr = pio_inq_varid(File, 'nextsw_cday', vardesc)
         765           0 :    ierr = pio_get_var(File, vardesc, temp_var)
         766           0 :    nextsw_cday = temp_var
         767             :
         768             : end subroutine radiation_read_restart
         769             :
         770             : !===============================================================================
         771             :
         772     5981472 : subroutine radiation_tend( &
         773             :    state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)
         774             :
         775             :    !-----------------------------------------------------------------------
         776             :    ! Driver for radiation computation.
         777             :    !
         778             :    ! NOTE: Radiation uses cgs units, so conversions must be done from
         779             :    !       model fields to radiation fields.
         780             :    !-----------------------------------------------------------------------
         781             :
         782             :    use physics_buffer,      only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
         783             :    use phys_grid,           only: get_rlat_all_p, get_rlon_all_p
         784             :    use physics_types,       only: physics_state, physics_ptend
         785             :    use time_manager,        only: get_curr_calday
         786             :    use radheat,             only: radheat_tend
         787             :    use physconst,           only: cpair, stebol
         788             :    use radconstants,        only: nlwbands, nswbands
         789             :    use radsw,               only: radcswmx
         790             :    use radlw,               only: radclwmx
         791             :    use rad_constituents,    only: rad_cnst_get_gas, rad_cnst_out
         792             :    use aer_rad_props,       only: aer_rad_props_sw, aer_rad_props_lw
         793             :    use interpolate_data,    only: vertinterp
         794             :    use radiation_data,      only: rad_data_write
         795             :    use cloud_cover_diags,   only: cloud_cover_diags_out
         796             :    use tropopause,          only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
         797             :    use orbit,               only: zenith
         798             :
         799             :    ! Arguments
         800             :    type(physics_state),       target, intent(in)    :: state
         801             :    type(physics_ptend),               intent(out)   :: ptend
         802             :    type(physics_buffer_desc),         pointer       :: pbuf(:)
         803             :    type(cam_out_t),                   intent(inout) :: cam_out
         804             :    type(cam_in_t),                    intent(in)    :: cam_in
         805             :    real(r8),                          intent(out)   :: net_flx(pcols)
         806             :    type(rad_out_t), target, optional, intent(out)   :: rd_out
         807             :
         808             :    ! Local variables
         809             :    type(rad_out_t), pointer :: rd  ! allow rd_out to be optional by allocating a local object
         810             :                                    ! if the argument is not present
         811             :
         812             :    integer :: i, k
         813             :    integer :: lchnk, ncol
         814             :
         815             :    logical :: dosw, dolw, doabsems
         816     1495368 :    integer, pointer :: nmxrgn(:)              ! pbuf pointer to Number of maximally overlapped regions
         817     1495368 :    real(r8),pointer :: pmxrgn(:,:)            ! Maximum values of pressure for each
         818             :                                               !    maximally overlapped region.
         819             :                                               !    0->pmxrgn(i,1) is range of pressure for
         820             :                                               !    1st region,pmxrgn(i,1)->pmxrgn(i,2) for
         821             :                                               !    2nd region, etc
         822             :
         823     1495368 :    real(r8),pointer :: emis(:,:)              ! Cloud longwave emissivity
         824     1495368 :    real(r8),pointer :: cldtau(:,:)            ! Cloud longwave optical depth
         825     1495368 :    real(r8),pointer :: cicewp(:,:)            ! in-cloud cloud ice water path
         826     1495368 :    real(r8),pointer :: cliqwp(:,:)            ! in-cloud cloud liquid water path
         827             :
         828             :    real(r8) :: cltot(pcols)                   ! Diagnostic total cloud cover
         829             :    real(r8) :: cllow(pcols)                   !       "     low  cloud cover
         830             :    real(r8) :: clmed(pcols)                   !       "     mid  cloud cover
         831             :    real(r8) :: clhgh(pcols)                   !       "     hgh  cloud cover
         832             :
         833             :    real(r8) :: ftem(pcols,pver)               ! Temporary workspace for outfld variables
         834             :
         835             :    integer :: itim_old
         836     1495368 :    real(r8), pointer, dimension(:,:) :: rel     ! liquid effective drop radius (microns)
         837     1495368 :    real(r8), pointer, dimension(:,:) :: rei     ! ice effective drop size (microns)
         838     1495368 :    real(r8), pointer, dimension(:,:) :: cld     ! cloud fraction
         839     1495368 :    real(r8), pointer, dimension(:,:) :: qrs     ! shortwave radiative heating rate
         840     1495368 :    real(r8), pointer, dimension(:,:) :: qrl     ! longwave  radiative heating rate
         841             :
         842             :    real(r8) :: calday                        ! current calendar day
         843             :    real(r8) :: clat(pcols)                   ! current latitudes(radians)
         844             :    real(r8) :: clon(pcols)                   ! current longitudes(radians)
         845             :    real(r8) :: coszrs(pcols)                 ! Cosine solar zenith angle
         846             :
         847             :    real(r8) :: fns(pcols,pverp)     ! net shortwave flux
         848             :    real(r8) :: fcns(pcols,pverp)    ! net clear-sky shortwave flux
         849             :    real(r8) :: fnl(pcols,pverp)     ! net longwave flux
         850             :    real(r8) :: fcnl(pcols,pverp)    ! net clear-sky longwave flux
         851             :
         852             :    ! This is used by the chemistry.
         853     1495368 :    real(r8), pointer :: fsds(:)  ! Surface solar down flux
         854             :
         855             :    ! This is used for the energy checker and the Eulerian dycore.
         856     1495368 :    real(r8), pointer :: fsns(:)  ! Surface solar absorbed flux
         857     1495368 :    real(r8), pointer :: fsnt(:)  ! Net column abs solar flux at model top
         858     1495368 :    real(r8), pointer :: flns(:)  ! Srf longwave cooling (up-down) flux
         859     1495368 :    real(r8), pointer :: flnt(:)  ! Net outgoing lw flux at model top
         860             :
         861             :    real(r8) :: pbr(pcols,pver)      ! Model mid-level pressures (dynes/cm2)
         862             :    real(r8) :: pnm(pcols,pverp)     ! Model interface pressures (dynes/cm2)
         863             :    real(r8) :: eccf                 ! Earth/sun distance factor
         864             :    real(r8) :: lwupcgs(pcols)       ! Upward longwave flux in cgs units
         865             :
         866     1495368 :    real(r8), pointer, dimension(:,:) :: n2o    ! nitrous oxide mass mixing ratio
         867     1495368 :    real(r8), pointer, dimension(:,:) :: ch4    ! methane mass mixing ratio
         868     1495368 :    real(r8), pointer, dimension(:,:) :: cfc11  ! cfc11 mass mixing ratio
         869     1495368 :    real(r8), pointer, dimension(:,:) :: cfc12  ! cfc12 mass mixing ratio
         870     1495368 :    real(r8), pointer, dimension(:,:) :: o3     ! Ozone mass mixing ratio
         871     1495368 :    real(r8), pointer, dimension(:,:) :: o2     ! Oxygen mass mixing ratio
         872             :    real(r8), dimension(pcols) :: o2_col        ! column oxygen mmr
         873     1495368 :    real(r8), pointer, dimension(:,:) :: co2    ! co2   mass mixing ratio
         874             :    real(r8), dimension(pcols) :: co2_col_mean  ! co2 column mean mmr
         875     1495368 :    real(r8), pointer, dimension(:,:) :: sp_hum ! specific humidity
         876             :
         877             :    ! Aerosol shortwave radiative properties
         878             :    real(r8) :: aer_tau    (pcols,0:pver,nswbands) ! aerosol extinction optical depth
         879             :    real(r8) :: aer_tau_w  (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
         880             :    real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * w * tau
         881             :    real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
         882             :
         883             :    ! Aerosol longwave absorption optical depth
         884             :    real(r8) :: odap_aer(pcols,pver,nlwbands)
         885             :
         886             :    ! Gathered indicies of day and night columns
         887             :    !  chunk_column_index = IdxDay(daylight_column_index)
         888             :    integer :: Nday                      ! Number of daylight columns
         889             :    integer :: Nnite                     ! Number of night columns
         890             :    integer, dimension(pcols) :: IdxDay  ! Indicies of daylight coumns
         891             :    integer, dimension(pcols) :: IdxNite ! Indicies of night coumns
         892             :
         893             :    character(*), parameter :: name = 'radiation_tend'
         894             :
         895             :    ! tropopause diagnostic
         896             :    integer :: troplev(pcols)
         897             :    real(r8):: p_trop(pcols)
         898             :
         899             :    logical :: write_output ! switch for outfld calls
         900             :    !----------------------------------------------------------------------
         901             :
         902     1495368 :    lchnk = state%lchnk
         903     1495368 :    ncol = state%ncol
         904             :
         905     1495368 :    calday = get_curr_calday()
         906             :
         907     1495368 :    if (present(rd_out)) then
         908             :       rd => rd_out
         909             :       write_output = .false.
         910             :    else
         911     1495368 :       allocate(rd)
         912             :       write_output=.true.
         913             :    end if
         914             :
         915     1495368 :    itim_old = pbuf_old_tim_idx()
         916     5981472 :    call pbuf_get_field(pbuf, cld_idx,    cld,    start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         917             :
         918     1495368 :    call pbuf_get_field(pbuf, qrs_idx,qrs)
         919     1495368 :    call pbuf_get_field(pbuf, qrl_idx,qrl)
         920             :
         921     1495368 :    call pbuf_get_field(pbuf, fsds_idx, fsds)
         922             :
         923     1495368 :    call pbuf_get_field(pbuf, fsns_idx, fsns)
         924     1495368 :    call pbuf_get_field(pbuf, fsnt_idx, fsnt)
         925     1495368 :    call pbuf_get_field(pbuf, flns_idx, flns)
         926     1495368 :    call pbuf_get_field(pbuf, flnt_idx, flnt)
         927             :
         928     1495368 :    call pbuf_get_field(pbuf, rel_idx, rel)
         929     1495368 :    call pbuf_get_field(pbuf, rei_idx, rei)
         930             :
         931             :    !  For CRM, make cloud equal to input observations:
         932     1495368 :    if (single_column.and.scm_crm_mode.and.have_cld) then
         933           0 :       do k = 1,pver
         934           0 :          cld(:ncol,k)= cldobs(k)
         935             :       enddo
         936             :    endif
         937             :
         938             :    ! Cosine solar zenith angle for current time step
         939     1495368 :    call get_rlat_all_p(lchnk, ncol, clat)
         940     1495368 :    call get_rlon_all_p(lchnk, ncol, clon)
         941     1495368 :    if (use_rad_uniform_angle) then
         942           0 :       call zenith (calday, clat, clon, coszrs, ncol, dt_avg, uniform_angle=rad_uniform_angle)
         943             :    else
         944     1495368 :       call zenith (calday, clat, clon, coszrs, ncol, dt_avg)
         945             :    end if
         946             :
         947             :    ! Gather night/day column indices.
         948     1495368 :    Nday = 0
         949     1495368 :    Nnite = 0
         950    24969168 :    do i = 1, ncol
         951    24969168 :       if ( coszrs(i) > 0.0_r8 ) then
         952    11736900 :          Nday = Nday + 1
         953    11736900 :          IdxDay(Nday) = i
         954             :       else
         955    11736900 :          Nnite = Nnite + 1
         956    11736900 :          IdxNite(Nnite) = i
         957             :       end if
         958             :    end do
         959             :
         960     1495368 :    dosw     = radiation_do('sw')      ! do shortwave heating calc this timestep?
         961     1495368 :    dolw     = radiation_do('lw')      ! do longwave heating calc this timestep?
         962             :
         963     1495368 :    doabsems = radiation_do('absems')  ! do absorptivity/emissivity calc this timestep?
         964             :
         965             :    ! Get time of next radiation calculation - albedos will need to be
         966             :    ! calculated by each surface model at this time
         967     1495368 :    nextsw_cday = radiation_nextsw_cday()
         968             :
         969     1495368 :    if (dosw .or. dolw) then
         970             :
         971             :       ! pbuf cloud properties set in cloud_diagnostics
         972      749232 :       call pbuf_get_field(pbuf, cicewp_idx, cicewp)
         973      749232 :       call pbuf_get_field(pbuf, cliqwp_idx, cliqwp)
         974      749232 :       call pbuf_get_field(pbuf, cldemis_idx, emis)
         975             :
         976      749232 :       call pbuf_get_field(pbuf, cldtau_idx, cldtau)
         977             :
         978      749232 :       call pbuf_get_field(pbuf, pmxrgn_idx, pmxrgn)
         979      749232 :       call pbuf_get_field(pbuf, nmxrgn_idx, nmxrgn)
         980             :
         981             :       ! For CRM, make cloud liquid water path equal to input observations
         982      749232 :       if(single_column.and.scm_crm_mode.and.have_clwp)then
         983           0 :          do k=1,pver
         984           0 :             cliqwp(:ncol,k) = clwpobs(k)
         985             :          end do
         986             :       endif
         987             :
         988             :       ! Get specific humidity
         989      749232 :       call rad_cnst_get_gas(0,'H2O', state, pbuf,  sp_hum)
         990             :
         991             :       ! Get ozone mass mixing ratio.
         992      749232 :       call rad_cnst_get_gas(0,'O3',  state, pbuf,  o3)
         993             :
         994             :       ! Get CO2 mass mixing ratio and compute column mean values
         995      749232 :       call rad_cnst_get_gas(0,'CO2', state, pbuf,  co2)
         996      749232 :       call calc_col_mean(state, co2, co2_col_mean)
         997             :
         998             :       ! construct cgs unit reps of pmid and pint and get "eccf" - earthsundistancefactor
         999      749232 :       call radinp(ncol, state%pmid, state%pint, pbr, pnm, eccf)
        1000             :
        1001             :       ! Solar radiation computation
        1002             :
        1003      749232 :       if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
        1004           0 :          call tropopause_find_cam(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
        1005             :       endif
        1006             :
        1007      749232 :       if (dosw) then
        1008             :
        1009      749232 :          call t_startf('rad_sw')
        1010             :
        1011             :          ! Get Oxygen mass mixing ratio.
        1012      749232 :          call rad_cnst_get_gas(0,'O2', state, pbuf,  o2)
        1013      749232 :          call calc_col_mean(state, o2, o2_col)
        1014             :
        1015             :          ! Get aerosol radiative properties.
        1016      749232 :          call t_startf('aero_optics_sw')
        1017             :          call aer_rad_props_sw(0, state, pbuf,  nnite, idxnite, &
        1018      749232 :             aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
        1019      749232 :          call t_stopf('aero_optics_sw')
        1020             :
        1021             :          call radcswmx(lchnk, &
        1022             :                ncol,       pnm,        pbr,        sp_hum,     o3,         &
        1023             :                o2_col,     cld,        cicewp,     cliqwp,     rel,        &
        1024             :                rei,        eccf,       coszrs,     rd%solin,      &
        1025             :                cam_in%asdir, cam_in%asdif, cam_in%aldir, cam_in%aldif, nmxrgn, &
        1026             :                pmxrgn,     qrs,        rd%qrsc,       fsnt,       rd%fsntc,      rd%fsdtoa, &
        1027             :                rd%fsntoa,     rd%fsutoa,     rd%fsntoac,    rd%fsnirt,     rd%fsnrtc,     rd%fsnirtsq,   &
        1028             :                fsns,       rd%fsnsc,      rd%fsdsc,      fsds,       cam_out%sols, &
        1029             :                cam_out%soll, cam_out%solsd, cam_out%solld, fns, fcns,      &
        1030             :                Nday,       Nnite,      IdxDay,     IdxNite,    co2_col_mean, &
        1031      749232 :                aer_tau,    aer_tau_w,  aer_tau_w_g, aer_tau_w_f , rd%liq_icld_vistau, rd%ice_icld_vistau  )
        1032             :
        1033      749232 :          call t_stopf('rad_sw')
        1034             :
        1035             :          !  Output net fluxes at 200 mb
        1036      749232 :          call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c)
        1037      749232 :          call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns,  rd%fsn200)
        1038      749232 :          if (hist_fld_active('FSNR')) then
        1039           0 :             do i = 1,ncol
        1040           0 :                call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
        1041           0 :                rd%fsnr(i) = rd%fsnr(i)*cgs2mks
        1042             :             enddo
        1043             :          else
        1044    12736944 :             rd%fsnr(:) = 0._r8
        1045             :          endif
        1046             :
        1047             :          ! Convert units of shortwave fields needed by rest of model from CGS to MKS
        1048             :
        1049    12510432 :          do i=1,ncol
        1050    11761200 :             rd%solin(i)   = rd%solin(i)   *cgs2mks
        1051    11761200 :             fsds(i)       = fsds(i)       *cgs2mks
        1052    11761200 :             rd%fsnirt(i)  = rd%fsnirt(i)  *cgs2mks
        1053    11761200 :             rd%fsnrtc(i)  = rd%fsnrtc(i)  *cgs2mks
        1054    11761200 :             rd%fsnirtsq(i)= rd%fsnirtsq(i)*cgs2mks
        1055    11761200 :             fsnt(i)       = fsnt(i)       *cgs2mks
        1056    11761200 :             rd%fsdtoa(i)  = rd%fsdtoa(i)  *cgs2mks
        1057    11761200 :             fsns(i)       = fsns(i)       *cgs2mks
        1058    11761200 :             rd%fsntc(i)   = rd%fsntc(i)   *cgs2mks
        1059    11761200 :             rd%fsnsc(i)   = rd%fsnsc(i)   *cgs2mks
        1060    11761200 :             rd%fsdsc(i)   = rd%fsdsc(i)   *cgs2mks
        1061    11761200 :             rd%fsntoa(i)  = rd%fsntoa(i)  *cgs2mks
        1062    11761200 :             rd%fsutoa(i)  = rd%fsutoa(i)  *cgs2mks
        1063    11761200 :             rd%fsntoac(i) = rd%fsntoac(i) *cgs2mks
        1064    11761200 :             rd%fsn200(i)  = rd%fsn200(i)  *cgs2mks
        1065    11761200 :             rd%fsn200c(i) = rd%fsn200c(i) *cgs2mks
        1066    12510432 :             rd%swcf(i)    = rd%fsntoa(i) - rd%fsntoac(i)
        1067             :          end do
        1068             :
        1069             :          ! initialize tau_cld_vistau and tau_icld_vistau as fillvalue, they will stay fillvalue for night columns
        1070   331909776 :          rd%tot_icld_vistau(1:pcols,1:pver) = fillvalue
        1071   331909776 :          rd%tot_cld_vistau(1:pcols,1:pver)  = fillvalue
        1072             :
        1073             :          ! only do calcs for tot_cld_vistau and tot_icld_vistau on daytime columns
        1074     6629832 :          do i=1,Nday
        1075             :             ! sum the water and ice optical depths to get total in-cloud cloud optical depth
        1076    11761200 :             rd%tot_icld_vistau(IdxDay(i),1:pver) = rd%liq_icld_vistau(IdxDay(i),1:pver) + &
        1077   170537400 :                                                    rd%ice_icld_vistau(IdxDay(i),1:pver)
        1078             :
        1079             :             ! sum wat and ice, multiply by cloud fraction to get grid-box value
        1080             :             rd%tot_cld_vistau(IdxDay(i),1:pver) = (rd%liq_icld_vistau(IdxDay(i),1:pver) + &
        1081   312421032 :                                                    rd%ice_icld_vistau(IdxDay(i),1:pver))*cld(IdxDay(i),1:pver)
        1082             :          end do
        1083             :
        1084             :          ! add fillvalue for night columns
        1085     6629832 :          do i = 1, Nnite
        1086   158776200 :             rd%liq_icld_vistau(IdxNite(i),:)  = fillvalue
        1087   159525432 :             rd%ice_icld_vistau(IdxNite(i),:)  = fillvalue
        1088             :          end do
        1089             :
        1090      749232 :          if (write_output) call radiation_output_sw(state, rd, cam_out, fsns, fsnt, fsds, qrs)
        1091             :
        1092             :       end if   ! dosw
        1093             :
        1094             :       ! Longwave radiation computation
        1095             :
        1096      749232 :       if (dolw) then
        1097             :
        1098      749232 :          call t_startf("rad_lw")
        1099             :
        1100             :          ! Convert upward longwave flux units to CGS
        1101             :
        1102    12510432 :          do i=1,ncol
        1103    11761200 :             lwupcgs(i) = cam_in%lwup(i)*1000._r8
        1104    11761200 :             if (single_column .and. scm_crm_mode .and. have_tg) &
        1105      749232 :                lwupcgs(i) = 1000*stebol*tground(1)**4
        1106             :          end do
        1107             :
        1108             :          ! Get gas phase constituents.
        1109      749232 :          call rad_cnst_get_gas(0,'N2O',   state, pbuf,  n2o)
        1110      749232 :          call rad_cnst_get_gas(0,'CH4',   state, pbuf,  ch4)
        1111      749232 :          call rad_cnst_get_gas(0,'CFC11', state, pbuf,  cfc11)
        1112      749232 :          call rad_cnst_get_gas(0,'CFC12', state, pbuf,  cfc12)
        1113             :
        1114             :          ! absems requires lw absorption optical depth and transmission through aerosols
        1115      749232 :          call t_startf('aero_optics_lw')
        1116      749232 :          if (doabsems) call aer_rad_props_lw(0, state, pbuf, odap_aer)
        1117      749232 :          call t_stopf('aero_optics_lw')
        1118             :
        1119             :          call radclwmx(lchnk, ncol, doabsems, &
        1120             :                  lwupcgs, state%t, sp_hum, o3, pbr, &
        1121             :                  pnm, state%lnpmid, state%lnpint, n2o, ch4, &
        1122             :                  cfc11, cfc12, cld, emis, pmxrgn, &
        1123             :                  nmxrgn, qrl, rd%qrlc, flns, flnt, rd%flnsc, &
        1124             :                  rd%flntc, cam_out%flwds, rd%fldsc, rd%flut, rd%flutc, &
        1125      749232 :                  fnl, fcnl, co2_col_mean, odap_aer)
        1126             :
        1127      749232 :          call t_stopf("rad_lw")
        1128             :
        1129             :          !  Output fluxes at 200 mb
        1130      749232 :          call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl,  rd%fln200)
        1131      749232 :          call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
        1132      749232 :          if (hist_fld_active('FLNR')) then
        1133           0 :             do i = 1,ncol
        1134           0 :                call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
        1135             :             enddo
        1136             :          else
        1137    12736944 :             rd%flnr(:) = 0._r8
        1138             :          endif
        1139             :
        1140             :          ! Convert units of longwave fields needed by rest of model from CGS to MKS
        1141             :
        1142    12510432 :          do i = 1, ncol
        1143    11761200 :             flnt(i)          = flnt(i)            *cgs2mks
        1144    11761200 :             rd%flut(i)       = rd%flut(i)         *cgs2mks
        1145    11761200 :             rd%flutc(i)      = rd%flutc(i)        *cgs2mks
        1146    11761200 :             rd%lwcf(i)       = rd%flutc(i) - rd%flut(i)
        1147    11761200 :             flns(i)          = flns(i)            *cgs2mks
        1148    11761200 :             rd%fldsc(i)      = rd%fldsc(i)        *cgs2mks
        1149    11761200 :             rd%flntc(i)      = rd%flntc(i)        *cgs2mks
        1150    11761200 :             rd%fln200(i)     = rd%fln200(i)       *cgs2mks
        1151    11761200 :             rd%fln200c(i)    = rd%fln200c(i)      *cgs2mks
        1152    11761200 :             rd%flnsc(i)      = rd%flnsc(i)        *cgs2mks
        1153    11761200 :             cam_out%flwds(i) = cam_out%flwds(i)   *cgs2mks
        1154    12510432 :             rd%flnr(i)       = rd%flnr(i)         *cgs2mks
        1155             :          end do
        1156             :
        1157      749232 :          if (write_output) call radiation_output_lw(state,  rd, cam_out, flns, flnt, qrl)
        1158             :
        1159             :       end if  ! dolw
        1160             :
        1161             :       ! Output aerosol mmr
        1162      749232 :       if (write_output) call rad_cnst_out(0, state, pbuf)
        1163             :
        1164             :       ! Cloud cover diagnostics
        1165             :       ! radsw can change pmxrgn and nmxrgn so cldsav needs to follow radsw
        1166      749232 :       if (write_output) call cloud_cover_diags_out(lchnk, ncol, cld, state%pmid, nmxrgn, pmxrgn )
        1167             :
        1168             :    else   !  if (dosw .or. dolw) then
        1169             :
        1170             :       ! convert radiative heating rates from Q*dp to Q for energy conservation
        1171    20145672 :       do k =1 , pver
        1172   324673272 :          do i = 1, ncol
        1173   304527600 :             qrs(i,k) = qrs(i,k)/state%pdel(i,k)
        1174   323927136 :             qrl(i,k) = qrl(i,k)/state%pdel(i,k)
        1175             :          end do
        1176             :       end do
        1177             :
        1178             :    end if   !  if (dosw .or. dolw) then
        1179             :
        1180             :    ! output rad inputs and resulting heating rates
        1181     1495368 :    call rad_data_write( pbuf, state, cam_in, coszrs )
        1182             :
        1183             :    ! Compute net radiative heating tendency
        1184             :    call radheat_tend(state, pbuf,  ptend, qrl, qrs, fsns, &
        1185     1495368 :                      fsnt, flns, flnt, cam_in%asdir, net_flx)
        1186             :
        1187     1495368 :    if (write_output) then
        1188             :       ! Compute heating rate for dtheta/dt
        1189    40374936 :       do k=1,pver
        1190   650693736 :          do i=1,ncol
        1191   649198368 :             ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
        1192             :          end do
        1193             :       end do
        1194     1495368 :       call outfld('HR      ',ftem    ,pcols   ,lchnk   )
        1195             :    end if
        1196             :
        1197             :    ! convert radiative heating rates to Q*dp for energy conservation
        1198    40374936 :    do k =1 , pver
        1199   650693736 :       do i = 1, ncol
        1200   610318800 :          qrs(i,k) = qrs(i,k)*state%pdel(i,k)
        1201   649198368 :          qrl(i,k) = qrl(i,k)*state%pdel(i,k)
        1202             :       end do
        1203             :    end do
        1204             :
        1205    24969168 :    cam_out%netsw(:ncol) = fsns(:ncol)
        1206             :
        1207     1495368 :    if (.not. present(rd_out)) then
        1208     1495368 :       deallocate(rd)
        1209             :    end if
        1210     2990736 : end subroutine radiation_tend
        1211             :
        1212             : !===============================================================================
        1213             :
        1214      749232 : subroutine radiation_output_sw(state, rd, cam_out, fsns, fsnt, fsds, qrs)
        1215             :
        1216             :    ! Dump shortwave radiation information to history buffer (diagnostics)
        1217             :
        1218             :    type(physics_state), intent(in) :: state
        1219             :    type(rad_out_t),     intent(in) :: rd
        1220             :    type(cam_out_t),     intent(in) :: cam_out
        1221             :    real(r8),            intent(in) :: fsns(pcols) ! Surface solar absorbed flux
        1222             :    real(r8),            intent(in) :: fsnt(pcols) ! Net column abs solar flux at model top
        1223             :    real(r8),            intent(in) :: fsds(pcols) ! Surface solar down flux
        1224             :    real(r8),            pointer    :: qrs(:,:)    ! shortwave radiative heating rate
        1225             :
        1226             :    ! Local variables
        1227             :    integer :: lchnk, ncol
        1228             :    real(r8) :: ftem(pcols,pver)               ! Temporary workspace for outfld variables
        1229             :    !----------------------------------------------------------------------
        1230             :
        1231      749232 :    lchnk = state%lchnk
        1232      749232 :    ncol = state%ncol
        1233             :
        1234   326020464 :    ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair
        1235      749232 :    call outfld('QRS     ',ftem  ,pcols,lchnk)
        1236             :
        1237   326020464 :    ftem(:ncol,:pver) = rd%qrsc(:ncol,:pver)/cpair
        1238      749232 :    call outfld('QRSC    ',ftem  ,pcols,lchnk)
        1239             :
        1240      749232 :    call outfld('SOLIN   ',rd%solin      ,pcols,lchnk)
        1241      749232 :    call outfld('FSDS    ',fsds       ,pcols,lchnk)
        1242      749232 :    call outfld('FSNIRTOA',rd%fsnirt     ,pcols,lchnk)
        1243      749232 :    call outfld('FSNRTOAC',rd%fsnrtc     ,pcols,lchnk)
        1244      749232 :    call outfld('FSNRTOAS',rd%fsnirtsq   ,pcols,lchnk)
        1245      749232 :    call outfld('FSNT    ',fsnt       ,pcols,lchnk)
        1246      749232 :    call outfld('FSDTOA  ',rd%fsdtoa     ,pcols,lchnk)
        1247      749232 :    call outfld('FSNS    ',fsns       ,pcols,lchnk)
        1248      749232 :    call outfld('FSNTC   ',rd%fsntc      ,pcols,lchnk)
        1249      749232 :    call outfld('FSNSC   ',rd%fsnsc      ,pcols,lchnk)
        1250      749232 :    call outfld('FSDSC   ',rd%fsdsc      ,pcols,lchnk)
        1251      749232 :    call outfld('FSNTOA  ',rd%fsntoa     ,pcols,lchnk)
        1252      749232 :    call outfld('FSUTOA  ',rd%fsutoa     ,pcols,lchnk)
        1253      749232 :    call outfld('FSNTOAC ',rd%fsntoac    ,pcols,lchnk)
        1254      749232 :    call outfld('SOLS    ',cam_out%sols  ,pcols,lchnk)
        1255      749232 :    call outfld('SOLL    ',cam_out%soll  ,pcols,lchnk)
        1256      749232 :    call outfld('SOLSD   ',cam_out%solsd ,pcols,lchnk)
        1257      749232 :    call outfld('SOLLD   ',cam_out%solld ,pcols,lchnk)
        1258      749232 :    call outfld('FSN200  ',rd%fsn200     ,pcols,lchnk)
        1259      749232 :    call outfld('FSN200C ',rd%fsn200c    ,pcols,lchnk)
        1260      749232 :    call outfld('FSNR'    ,rd%fsnr       ,pcols,lchnk)
        1261      749232 :    call outfld('SWCF    ',rd%swcf       ,pcols,lchnk)
        1262             :
        1263      749232 :    call outfld('TOT_CLD_VISTAU    ',rd%tot_cld_vistau     ,pcols,lchnk)
        1264      749232 :    call outfld('TOT_ICLD_VISTAU   ',rd%tot_icld_vistau    ,pcols,lchnk)
        1265      749232 :    call outfld('LIQ_ICLD_VISTAU   ',rd%liq_icld_vistau ,pcols,lchnk)
        1266      749232 :    call outfld('ICE_ICLD_VISTAU   ',rd%ice_icld_vistau ,pcols,lchnk)
        1267             :
        1268     1495368 : end subroutine radiation_output_sw
        1269             :
        1270             : !===============================================================================
        1271             :
        1272      749232 : subroutine radiation_output_lw(state,  rd, cam_out, flns, flnt, qrl)
        1273             :
        1274             :    ! Dump longwave radiation information to history tape buffer (diagnostics)
        1275             :
        1276             :    type(physics_state), intent(in) :: state
        1277             :    type(rad_out_t),     intent(in) :: rd
        1278             :    type(cam_out_t),     intent(in) :: cam_out
        1279             :    real(r8),            intent(in) :: flns(pcols) ! Srf longwave cooling (up-down) flux
        1280             :    real(r8),            intent(in) :: flnt(pcols) ! Net outgoing lw flux at model top
        1281             :    real(r8),            pointer    :: qrl(:,:)    ! longwave  radiative heating rate
        1282             :
        1283             :    ! Local variables
        1284             :    integer :: lchnk, ncol
        1285             :    real(r8) :: ftem(pcols,pver)
        1286             :    !----------------------------------------------------------------------
        1287             :
        1288      749232 :    lchnk = state%lchnk
        1289      749232 :    ncol = state%ncol
        1290             :
        1291   326020464 :    call outfld('QRL     ',qrl(:ncol,:)/cpair,ncol,lchnk)
        1292   326020464 :    call outfld('QRLC    ',rd%qrlc(:ncol,:)/cpair,ncol,lchnk)
        1293      749232 :    call outfld('FLNT    ',flnt  ,pcols,lchnk)
        1294      749232 :    call outfld('FLUT    ',rd%flut  ,pcols,lchnk)
        1295      749232 :    call outfld('FLUTC   ',rd%flutc ,pcols,lchnk)
        1296      749232 :    call outfld('FLNTC   ',rd%flntc ,pcols,lchnk)
        1297      749232 :    call outfld('FLNS    ',flns  ,pcols,lchnk)
        1298      749232 :    call outfld('FLDS    ',cam_out%flwds ,pcols,lchnk)
        1299      749232 :    call outfld('FLNSC   ',rd%flnsc ,pcols,lchnk)
        1300      749232 :    call outfld('FLDSC   ',rd%fldsc ,pcols,lchnk)
        1301      749232 :    call outfld('LWCF    ',rd%lwcf  ,pcols,lchnk)
        1302      749232 :    call outfld('FLN200  ',rd%fln200,pcols,lchnk)
        1303      749232 :    call outfld('FLN200C ',rd%fln200c,pcols,lchnk)
        1304      749232 :    call outfld('FLNR '   ,rd%flnr,pcols,lchnk)
        1305             :
        1306      749232 : end subroutine radiation_output_lw
        1307             :
        1308             : !===============================================================================
        1309             :
        1310      749232 : subroutine radinp(ncol, pmid, pint, pmidrd, pintrd, eccf)
        1311             :
        1312             :    use shr_orb_mod
        1313             :    use time_manager, only: get_curr_calday
        1314             :
        1315             :    !------------------------------Arguments--------------------------------
        1316             :    integer, intent(in)   :: ncol                ! number of atmospheric columns
        1317             :
        1318             :    real(r8), intent(in)  :: pmid(pcols,pver)    ! Pressure at model mid-levels (pascals)
        1319             :    real(r8), intent(in)  :: pint(pcols,pverp)   ! Pressure at model interfaces (pascals)
        1320             :
        1321             :    real(r8), intent(out) :: pmidrd(pcols,pver)  ! Pressure at mid-levels (dynes/cm*2)
        1322             :    real(r8), intent(out) :: pintrd(pcols,pverp) ! Pressure at interfaces (dynes/cm*2)
        1323             :    real(r8), intent(out) :: eccf                ! Earth-sun distance factor
        1324             :
        1325             :    !---------------------------Local variables-----------------------------
        1326             :    integer :: i, k
        1327             :    real(r8) :: calday       ! current calendar day
        1328             :    real(r8) :: delta        ! Solar declination angle
        1329             :    !-----------------------------------------------------------------------
        1330             :
        1331      749232 :    calday = get_curr_calday()
        1332             :    call shr_orb_decl (calday  ,eccen     ,mvelpp  ,lambm0  ,obliqr  , &
        1333      749232 :                       delta   ,eccf)
        1334             :
        1335             :    ! Convert pressure from pascals to dynes/cm2
        1336    20229264 :    do k=1,pver
        1337   326020464 :       do i=1,ncol
        1338   305791200 :          pmidrd(i,k) = pmid(i,k)*10.0_r8
        1339   325271232 :          pintrd(i,k) = pint(i,k)*10.0_r8
        1340             :       end do
        1341             :    end do
        1342    12510432 :    do i=1,ncol
        1343    12510432 :       pintrd(i,pverp) = pint(i,pverp)*10.0_r8
        1344             :    end do
        1345             :
        1346      749232 : end subroutine radinp
        1347             :
        1348             : !===============================================================================
        1349             :
        1350     1498464 : subroutine calc_col_mean(state, mmr_pointer, mean_value)
        1351             :
        1352             :    ! Compute the column mean.
        1353             :
        1354      749232 :    use cam_logfile,  only: iulog
        1355             :
        1356             :    type(physics_state),        intent(in)  :: state
        1357             :    real(r8), dimension(:,:),   pointer     :: mmr_pointer  ! mass mixing ratio (lev)
        1358             :    real(r8), dimension(pcols), intent(out) :: mean_value   ! column mean mmr
        1359             :
        1360             :    integer  :: i, k, ncol
        1361             :    real(r8) :: ptot(pcols)
        1362             :    !-----------------------------------------------------------------------
        1363             :
        1364     1498464 :    ncol         = state%ncol
        1365     1498464 :    mean_value   = 0.0_r8
        1366     1498464 :    ptot         = 0.0_r8
        1367             :
        1368    40458528 :    do k=1,pver
        1369   652040928 :       do i=1,ncol
        1370   611582400 :          mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k)
        1371   650542464 :          ptot(i)         = ptot(i) + state%pdeldry(i,k)
        1372             :       end do
        1373             :    end do
        1374    25020864 :    do i=1,ncol
        1375    25020864 :       mean_value(i) = mean_value(i) / ptot(i)
        1376             :    end do
        1377             :
        1378     1498464 : end subroutine calc_col_mean
        1379             :
        1380             : !===============================================================================
        1381             :
        1382           0 : end module radiation
        1383             :

Generated by: LCOV version 1.14