LCOV - code coverage report
Current view: | top level - physics/cam - aer_rad_props.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 157 | 245 | 64.1 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 6 | 8 | 75.0 %

              Line data    Source code
           1             : module aer_rad_props
           2             :
           3             : !------------------------------------------------------------------------------------------------
           4             : ! Converts aerosol masses to bulk optical properties for sw and lw radiation
           5             : ! computations.
           6             : !------------------------------------------------------------------------------------------------
           7             :
           8             : use shr_kind_mod,     only: r8 => shr_kind_r8
           9             : use ppgrid,           only: pcols, pver
          10             : use physconst,        only: rga
          11             : use physics_types,    only: physics_state
          12             :
          13             : use physics_buffer,   only: physics_buffer_desc
          14             : use radconstants,     only: nswbands, nlwbands, idx_sw_diag
          15             : use phys_prop,        only: nrh, ot_length
          16             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, &
          17             :                             rad_cnst_get_aer_props
          18             : use wv_saturation,    only: qsat
          19             : use aerosol_optics_cam,only: aerosol_optics_cam_init, aerosol_optics_cam_sw, aerosol_optics_cam_lw
          20             : use cam_history,      only: fieldname_len, addfld, outfld, add_default, horiz_only
          21             : use cam_history_support, only : fillvalue
          22             : ! Placed here due to PGI bug.
          23             : use ref_pres,         only: clim_modal_aero_top_lev
          24             :
          25             : use cam_abortutils,   only: endrun
          26             :
          27             : implicit none
          28             : private
          29             : save
          30             :
          31             : integer :: top_lev = 1
          32             :
          33             : public :: &
          34             :    aer_rad_props_init,        &
          35             :    aer_rad_props_sw,          & ! return SW optical props of aerosols
          36             :    aer_rad_props_lw             ! return LW optical props of aerosols
          37             :
          38             : ! Private data
          39             : character(len=fieldname_len), pointer :: odv_names(:)  ! outfld names for visible OD
          40             :
          41             :
          42             : !==============================================================================
          43             : contains
          44             : !==============================================================================
          45             :
          46        1536 : subroutine aer_rad_props_init()
          47             :    use phys_control, only: phys_getopts
          48             :
          49             :
          50             :    integer                    :: i
          51             :    integer                    :: numaerosols  ! number of aerosols
          52        1536 :    character(len=64), pointer :: aernames(:)  ! aerosol names
          53             :    logical                    :: history_amwg         ! output the variables used by the AMWG diag package
          54             :    logical                    :: history_aero_optics  ! Output aerosol optics diagnostics
          55             :    logical                    :: history_dust         ! Output dust diagnostics
          56             :    logical                    :: prog_modal_aero      ! Prognostic modal aerosols present
          57             :    integer                    :: nmodes          ! number of aerosol modes
          58             :
          59             :    !----------------------------------------------------------------------------
          60             :
          61             :    call phys_getopts( history_aero_optics_out    = history_aero_optics, &
          62             :                       history_amwg_out           = history_amwg,    &
          63             :                       history_dust_out           = history_dust,    &
          64        1536 :                       prog_modal_aero_out        = prog_modal_aero )
          65             :
          66             :    ! Limit modal aerosols with top_lev here.
          67        1536 :    if (prog_modal_aero) top_lev = clim_modal_aero_top_lev
          68             :
          69             :    call addfld ('AEROD_v ', horiz_only, 'A', '1', &
          70        1536 :       'Total Aerosol Optical Depth in visible band', flag_xyfill=.true.)
          71             :
          72             :    call addfld ('AODvstrt', horiz_only, 'A', '1', &
          73        1536 :       'Stratospheric Aerosol Optical Depth in visible band', flag_xyfill=.true.)
          74             :
          75             :    ! Contributions to AEROD_v from individual aerosols (climate species).
          76             :
          77             :    ! number of bulk aerosols in climate list
          78        1536 :    call rad_cnst_get_info(0, naero=numaerosols)
          79             :
          80             :    ! get names of bulk aerosols
          81        4608 :    allocate(aernames(numaerosols))
          82        1536 :    call rad_cnst_get_info(0, aernames=aernames, nmodes=nmodes)
          83             :
          84             :    ! diagnostic output for bulk aerosols
          85             :    ! create outfld names for visible OD
          86        4608 :    allocate(odv_names(numaerosols))
          87       21504 :    do i = 1, numaerosols
          88       19968 :       odv_names(i) = 'ODV_'//trim(aernames(i))
          89           0 :       call addfld (odv_names(i), horiz_only, 'A', '1', &
          90       21504 :          trim(aernames(i))//' optical depth in visible band', flag_xyfill=.true.)
          91             :    end do
          92             :
          93             :    ! Determine default fields
          94        1536 :    if (history_amwg .or. history_dust ) then
          95        1536 :       call add_default ('AEROD_v', 1, ' ')
          96             :    endif
          97             :
          98        1536 :    if ( history_aero_optics ) then
          99           0 :       call add_default ('AEROD_v', 1, ' ')
         100           0 :       do i = 1, numaerosols
         101           0 :          odv_names(i) = 'ODV_'//trim(aernames(i))
         102           0 :          call add_default (odv_names(i), 1, ' ')
         103             :       end do
         104             :     endif
         105             :
         106        1536 :     if (nmodes > 0) then
         107           0 :        call aerosol_optics_cam_init()
         108             :     end if
         109             :
         110        1536 :    deallocate(aernames)
         111             :
         112        1536 : end subroutine aer_rad_props_init
         113             :
         114             : !==============================================================================
         115             :
         116      746136 : subroutine aer_rad_props_sw(list_idx, state, pbuf,  nnite, idxnite, &
         117             :                             tau, tau_w, tau_w_g, tau_w_f)
         118             :
         119             :    ! Return bulk layer tau, omega, g, f for all spectral intervals.
         120             :
         121             :    use physics_buffer, only : physics_buffer_desc
         122             :    use tropopause,     only : tropopause_find_cam
         123             :    ! Arguments
         124             :    integer,             intent(in) :: list_idx      ! index of the climate or a diagnostic list
         125             :    type(physics_state), intent(in), target :: state
         126             :
         127             :    type(physics_buffer_desc), pointer :: pbuf(:)
         128             :    integer,             intent(in) :: nnite                ! number of night columns
         129             :    integer,             intent(in) :: idxnite(:)           ! local column indices of night columns
         130             :
         131             :    real(r8), intent(out) :: tau    (pcols,0:pver,nswbands) ! aerosol extinction optical depth
         132             :    real(r8), intent(out) :: tau_w  (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
         133             :    real(r8), intent(out) :: tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * tau * w
         134             :    real(r8), intent(out) :: tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * tau * w
         135             :
         136             :    ! Local variables
         137             :
         138             :    integer :: ncol
         139             :    integer :: lchnk
         140             :    integer :: k       ! index
         141             :    integer :: troplev(pcols)
         142             :
         143             :    ! optical props for each aerosol
         144             :    ! hygroscopic
         145      746136 :    real(r8), pointer :: h_ext(:,:)
         146      746136 :    real(r8), pointer :: h_ssa(:,:)
         147      746136 :    real(r8), pointer :: h_asm(:,:)
         148             :    ! non-hygroscopic
         149      746136 :    real(r8), pointer :: n_ext(:)
         150      746136 :    real(r8), pointer :: n_ssa(:)
         151      746136 :    real(r8), pointer :: n_asm(:)
         152      746136 :    real(r8), pointer :: n_scat(:)
         153      746136 :    real(r8), pointer :: n_ascat(:)
         154             :    ! radius-dependent
         155      746136 :    real(r8), pointer :: r_ext(:,:)    ! radius-dependent mass-specific extinction
         156      746136 :    real(r8), pointer :: r_scat(:,:)
         157      746136 :    real(r8), pointer :: r_ascat(:,:)
         158      746136 :    real(r8), pointer :: r_mu(:)       ! log(radius) domain variable for r_ext, r_scat, r_ascat
         159             :
         160             :    ! radiative properties for each aerosol
         161             :    real(r8) :: ta (pcols,pver,nswbands)
         162             :    real(r8) :: tw (pcols,pver,nswbands)
         163             :    real(r8) :: twf(pcols,pver,nswbands)
         164             :    real(r8) :: twg(pcols,pver,nswbands)
         165             :
         166             :    ! aerosol masses
         167      746136 :    real(r8), pointer :: aermmr(:,:)    ! mass mixing ratio of aerosols
         168             :    real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass
         169             :    real(r8) :: aermass(pcols,pver)     ! mass of aerosols
         170             :
         171             :    ! for table lookup into rh grid
         172             :    real(r8) :: es(pcols,pver)     ! saturation vapor pressure
         173             :    real(r8) :: qs(pcols,pver)     ! saturation specific humidity
         174             :    real(r8) :: rh(pcols,pver)
         175             :    real(r8) :: rhtrunc(pcols,pver)
         176             :    real(r8) :: wrh(pcols,pver)
         177             :    integer  :: krh(pcols,pver)
         178             :
         179             :    integer  :: numaerosols     ! number of bulk aerosols in climate/diagnostic list
         180             :    integer  :: nmodes          ! number of aerosol modes in climate/diagnostic list
         181             :    integer  :: iaerosol        ! index into bulk aerosol list
         182             :
         183             :    character(len=ot_length) :: opticstype       ! hygro or nonhygro
         184             :    character(len=16) :: pbuf_fld
         185             :    !-----------------------------------------------------------------------------
         186             :
         187      746136 :    ncol  = state%ncol
         188      746136 :    lchnk = state%lchnk
         189             :
         190             :    ! compute mixing ratio to mass conversion
         191    70136784 :    do k = 1, pver
         192  1159408584 :       mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k)
         193             :    enddo
         194             :
         195             :    ! initialize to conditions that would cause failure
         196 16703746632 :    tau     (:,:,:) = -100._r8
         197 16703746632 :    tau_w   (:,:,:) = -100._r8
         198 16703746632 :    tau_w_g (:,:,:) = -100._r8
         199 16703746632 :    tau_w_f (:,:,:) = -100._r8
         200             :
         201             :    ! top layer (ilev = 0) has no aerosol (ie tau = 0)
         202             :    ! also initialize rest of layers to accumulate od's
         203 16406888616 :    tau    (1:ncol,:,:) = 0._r8
         204 16406888616 :    tau_w  (1:ncol,:,:) = 0._r8
         205 16406888616 :    tau_w_g(1:ncol,:,:) = 0._r8
         206 16406888616 :    tau_w_f(1:ncol,:,:) = 0._r8
         207             :
         208             :    ! calculate relative humidity for table lookup into rh grid
         209    70136784 :    do k = 1, pver
         210    70136784 :       call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
         211             :    end do
         212  1159408584 :    rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver)
         213             :
         214  1159408584 :    rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8)
         215  1159408584 :    krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh
         216  1159408584 :    wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver)       ! (-) weighting on left side values
         217             :
         218             :    ! get number of bulk aerosols and number of modes in current list
         219      746136 :    call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes)
         220             :
         221             :    ! Contributions from modal aerosols.
         222      746136 :    if (nmodes > 0) then
         223             :       call aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, &
         224           0 :                                  tau, tau_w, tau_w_g, tau_w_f)
         225             :    else
         226 16406888616 :       tau    (1:ncol,:,:) = 0._r8
         227 16406888616 :       tau_w  (1:ncol,:,:) = 0._r8
         228 16406888616 :       tau_w_g(1:ncol,:,:) = 0._r8
         229 16406888616 :       tau_w_f(1:ncol,:,:) = 0._r8
         230             :    end if
         231             :
         232             :    !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         233      746136 :    troplev = 0
         234             :    !REMOVECAM_END
         235      746136 :    call tropopause_find_cam(state, troplev)
         236             :
         237             :    ! Contributions from bulk aerosols.
         238    10445904 :    do iaerosol = 1, numaerosols
         239             :
         240             :       ! get bulk aerosol mass mixing ratio
         241     9699768 :       call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf, aermmr)
         242     9699768 :       aermass(1:ncol,1:top_lev-1) = 0._r8
         243 15072311592 :       aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver)
         244             :
         245             :       ! get optics type
         246     9699768 :       call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype)
         247             :
         248    19399536 :       select case (trim(opticstype))
         249             :       case('hygro','hygroscopic','hygroscopi')
         250             :          ! get optical properties for hygroscopic aerosols
         251     4476816 :          call rad_cnst_get_aer_props(list_idx, iaerosol, sw_hygro_ext=h_ext, sw_hygro_ssa=h_ssa, sw_hygro_asm=h_asm)
         252     4476816 :          call get_hygro_rad_props(ncol, krh, wrh, aermass, h_ext, h_ssa, h_asm, ta, tw, twg, twf)
         253 97394797872 :          tau    (1:ncol,1:pver,:) = tau    (1:ncol,1:pver,:) + ta (1:ncol,:,:)
         254 97394797872 :          tau_w  (1:ncol,1:pver,:) = tau_w  (1:ncol,1:pver,:) + tw (1:ncol,:,:)
         255 97394797872 :          tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
         256 97394797872 :          tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
         257             :
         258             :       case('nonhygro','insoluble ')
         259             :          ! get optical properties for non-hygroscopic aerosols
         260             :          call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_ssa=n_ssa, &
         261     5222952 :                                      sw_nonhygro_asm=n_asm)
         262             :
         263     5222952 :          call get_nonhygro_rad_props(ncol, aermass, n_ext, n_ssa, n_asm, ta, tw, twg, twf)
         264 >11362*10^7 :          tau    (1:ncol,1:pver,:) = tau    (1:ncol,1:pver,:) + ta (1:ncol,:,:)
         265 >11362*10^7 :          tau_w  (1:ncol,1:pver,:) = tau_w  (1:ncol,1:pver,:) + tw (1:ncol,:,:)
         266 >11362*10^7 :          tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
         267 >11362*10^7 :          tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
         268             :
         269             :       case('volcanic')
         270             :          ! get optical properties for volcanic aerosols
         271             :          call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_scat=n_scat, &
         272           0 :                                      sw_nonhygro_ascat=n_ascat)
         273             :
         274           0 :          call get_volcanic_rad_props(ncol, aermass, n_ext, n_scat, n_ascat, ta, tw, twg, twf)
         275           0 :          tau    (1:ncol,1:pver,:) = tau    (1:ncol,1:pver,:) + ta (1:ncol,:,:)
         276           0 :          tau_w  (1:ncol,1:pver,:) = tau_w  (1:ncol,1:pver,:) + tw (1:ncol,:,:)
         277           0 :          tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
         278           0 :          tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
         279             :
         280             :       case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
         281           0 :          pbuf_fld = 'VOLC_RAD_GEOM '
         282           0 :          if (len_trim(opticstype)>15) then
         283           0 :             pbuf_fld = trim(pbuf_fld)//opticstype(16:16)
         284             :          endif
         285             :          ! get optical properties for volcanic aerosols
         286           0 :          call rad_cnst_get_aer_props(list_idx, iaerosol, r_sw_ext=r_ext, r_sw_scat=r_scat, r_sw_ascat=r_ascat, mu=r_mu)
         287           0 :          call get_volcanic_radius_rad_props(ncol, aermass, pbuf_fld, pbuf, r_ext, r_scat, r_ascat, r_mu, ta, tw, twg, twf)
         288           0 :          tau    (1:ncol,1:pver,:) = tau    (1:ncol,1:pver,:) + ta (1:ncol,:,:)
         289           0 :          tau_w  (1:ncol,1:pver,:) = tau_w  (1:ncol,1:pver,:) + tw (1:ncol,:,:)
         290           0 :          tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
         291           0 :          tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
         292             :
         293             :       case('zero')
         294             :          ! no effect of "zero" aerosols, so update nothing
         295             :       case default
         296    19399536 :          call endrun('aer_rad_props_sw: unsupported opticstype :'//trim(opticstype)//':')
         297             :       end select
         298             :
         299             :       ! diagnostic output of individual aerosol optical properties
         300             :       ! currently implemented for climate list only
         301    10445904 :       call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaerosol, ta(:,:,idx_sw_diag), list_idx, troplev)
         302             :
         303             :    enddo
         304             :
         305             :    ! diagnostic output of total aerosol optical properties
         306             :    ! currently implemented for climate list only
         307      746136 :    call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, 0, tau(:,:,idx_sw_diag), list_idx, troplev)
         308             :
         309     1492272 : end subroutine aer_rad_props_sw
         310             :
         311             : !==============================================================================
         312             :
         313      746136 : subroutine aer_rad_props_lw(list_idx, state, pbuf,  odap_aer)
         314             :
         315             :    ! Purpose: Compute aerosol transmissions needed in absorptivity/
         316             :    !    emissivity calculations
         317             :
         318             :    ! lw extinction is the same representation for all
         319             :    ! species.  If this changes, this routine will need to do something
         320             :    ! similar to the sw with routines like get_hygro_lw_abs
         321             :
         322      746136 :    use physics_buffer, only : pbuf_get_field, pbuf_get_index, physics_buffer_desc
         323             :
         324             :    ! Arguments
         325             :    integer,             intent(in)  :: list_idx                      ! index of the climate or a diagnostic list
         326             :    type(physics_state), intent(in), target :: state
         327             :
         328             :    type(physics_buffer_desc), pointer :: pbuf(:)
         329             :    real(r8),            intent(out) :: odap_aer(pcols,pver,nlwbands) ! [fraction] absorption optical depth, per layer
         330             :
         331             :    ! Local variables
         332             :
         333             :    integer :: bnd_idx     ! LW band index
         334             :    integer :: i           ! column index
         335             :    integer :: k           ! lev index
         336             :    integer :: ncol        ! number of columns
         337             :    integer :: numaerosols ! number of bulk aerosols in climate/diagnostic list
         338             :    integer :: nmodes      ! number of aerosol modes in climate/diagnostic list
         339             :    integer :: iaerosol    ! index into bulk aerosol list
         340             :    character(len=ot_length) :: opticstype       ! hygro or nonhygro
         341             :
         342             :    ! optical props for each aerosol
         343      746136 :    real(r8), pointer :: lw_abs(:)
         344      746136 :    real(r8), pointer :: lw_hygro_abs(:,:)
         345      746136 :    real(r8), pointer :: geometric_radius(:,:)
         346             :
         347             :    ! volcanic lookup table
         348      746136 :    real(r8), pointer :: r_lw_abs(:,:)  ! radius dependent mass-specific absorption coefficient
         349      746136 :    real(r8), pointer :: r_mu(:)        ! log(geometric_mean_radius) domain samples of r_lw_abs(:,:)
         350             :    integer  :: idx                     ! index to pbuf for geometric radius
         351             :    real(r8) :: mu(pcols,pver)          ! log(geometric_radius)
         352             :    real(r8) :: r_mu_min, r_mu_max, wmu, mutrunc
         353             :    integer  :: nmu, kmu
         354             :
         355             :    ! for table lookup into rh grid
         356             :    real(r8) :: es(pcols,pver)     ! saturation vapor pressure
         357             :    real(r8) :: qs(pcols,pver)     ! saturation specific humidity
         358             :    real(r8) :: rh(pcols,pver)
         359             :    real(r8) :: rhtrunc(pcols,pver)
         360             :    real(r8) :: wrh(pcols,pver)
         361             :    integer  :: krh(pcols,pver)
         362             :
         363             :    ! aerosol (vertical) mass path and extinction
         364             :    ! aerosol masses
         365      746136 :    real(r8), pointer :: aermmr(:,:)    ! mass mixing ratio of aerosols
         366             :    real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass
         367             :    real(r8) :: aermass(pcols,pver)     ! mass of aerosols
         368             :
         369             :    character(len=16) :: pbuf_fld
         370             :    !-----------------------------------------------------------------------------
         371             :
         372      746136 :    ncol = state%ncol
         373             :
         374             :    ! get number of bulk aerosols and number of modes in current list
         375      746136 :    call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes)
         376             :
         377             :    ! Contributions from modal aerosols.
         378      746136 :    if (nmodes > 0) then
         379           0 :       call aerosol_optics_cam_lw(list_idx, state, pbuf, odap_aer)
         380             :    else
         381      746136 :       odap_aer = 0._r8
         382             :    end if
         383             :
         384             :    ! Contributions from bulk aerosols.
         385      746136 :    if (numaerosols > 0) then
         386             :
         387             :       ! compute mixing ratio to mass conversion
         388    70136784 :       do k = 1, pver
         389  1159408584 :          mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k)
         390             :       end do
         391             :
         392             :       ! calculate relative humidity for table lookup into rh grid
         393    70136784 :       do k = 1, pver
         394    70136784 :          call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
         395             :       end do
         396  1159408584 :       rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver)
         397             :
         398  1159408584 :       rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8)
         399  1159408584 :       krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh
         400  1159408584 :       wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver)       ! (-) weighting on left side values
         401             :
         402             :    end if
         403             :
         404             :    ! Loop over bulk aerosols in list.
         405    10445904 :    do iaerosol = 1, numaerosols
         406             :
         407             :       ! get aerosol mass mixing ratio
         408     9699768 :       call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf,  aermmr)
         409     9699768 :       aermass(1:ncol,1:top_lev-1) = 0._r8
         410 15072311592 :       aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver)
         411             :
         412             :       ! get optics type
         413     9699768 :       call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype)
         414    20145672 :       select case (trim(opticstype))
         415             :       case('hygroscopic')
         416             :           ! get optical properties for hygroscopic aerosols
         417     4476816 :          call rad_cnst_get_aer_props(list_idx, iaerosol, lw_hygro_ext=lw_hygro_abs)
         418    76105872 :          do bnd_idx = 1, nlwbands
         419  6737608080 :             do k = 1, pver
         420 >11130*10^7 :                do i = 1, ncol
         421 >31371*10^7 :                   odap_aer(i, k, bnd_idx) = odap_aer(i, k, bnd_idx) + &
         422             :                        aermass(i, k) * &
         423           0 :                        ((1 + wrh(i,k)) * lw_hygro_abs(krh(i,k)+1,bnd_idx) &
         424 >42494*10^7 :                        - wrh(i,k)  * lw_hygro_abs(krh(i,k),  bnd_idx))
         425             :                end do
         426             :             end do
         427             :          end do
         428             :       case('insoluble','nonhygro','hygro','volcanic')
         429             :           ! get optical properties for hygroscopic aerosols
         430     5222952 :          call rad_cnst_get_aer_props(list_idx, iaerosol, lw_ext=lw_abs)
         431    88790184 :          do bnd_idx = 1, nlwbands
         432  7860542760 :             do k = 1, pver
         433 >12985*10^7 :                do i = 1, ncol
         434 >12977*10^7 :                   odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + lw_abs(bnd_idx)*aermass(i,k)
         435             :                end do
         436             :             end do
         437             :          end do
         438             :
         439             :       case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
         440           0 :          pbuf_fld = 'VOLC_RAD_GEOM '
         441           0 :          if (len_trim(opticstype)>15) then
         442           0 :             pbuf_fld = trim(pbuf_fld)//opticstype(16:16)
         443             :          endif
         444             :
         445             :          ! get optical properties for hygroscopic aerosols
         446           0 :          call rad_cnst_get_aer_props(list_idx, iaerosol, r_lw_abs=r_lw_abs, mu=r_mu)
         447             :          ! get microphysical properties for volcanic aerosols
         448           0 :          idx = pbuf_get_index(pbuf_fld)
         449           0 :          call pbuf_get_field(pbuf, idx, geometric_radius )
         450             :
         451             :          ! interpolate in radius
         452             :          ! caution: clip the table with no warning when outside bounds
         453           0 :          nmu = size(r_mu)
         454           0 :          r_mu_max = r_mu(nmu)
         455           0 :          r_mu_min = r_mu(1)
         456           0 :          do i = 1, ncol
         457           0 :             do k = 1, pver
         458           0 :                if(geometric_radius(i,k) > 0._r8) then
         459           0 :                   mu(i,k) = log(geometric_radius(i,k))
         460             :                else
         461           0 :                   mu(i,k) = 0._r8
         462             :                endif
         463           0 :                mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min)
         464           0 :                kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8)
         465           0 :                wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8)
         466           0 :                do bnd_idx = 1, nlwbands
         467           0 :                   odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + &
         468             :                      aermass(i,k) * &
         469           0 :                      ((1._r8 - wmu) * r_lw_abs(bnd_idx, kmu  ) + &
         470           0 :                      (wmu) * r_lw_abs(bnd_idx, kmu+1))
         471             :                end do
         472             :             end do
         473             :          end do
         474             :
         475             :        case('zero')
         476             :           ! zero aerosols types have no optical effect, so do nothing.
         477             :        case default
         478    19399536 :           call endrun('aer_rad_props_lw: unsupported opticstype: '//trim(opticstype))
         479             :        end select
         480             :     end do
         481             :
         482     1492272 : end subroutine aer_rad_props_lw
         483             :
         484             : !==============================================================================
         485             : ! Private methods
         486             : !==============================================================================
         487             :
         488     4476816 : subroutine get_hygro_rad_props(ncol, krh, wrh, mass, ext, ssa, asm, &
         489             :                                tau, tau_w, tau_w_g, tau_w_f)
         490             :
         491             :    ! Arguments
         492             :    integer,  intent(in) :: ncol
         493             :    integer,  intent(in) :: krh(pcols,pver)  ! index for linear interpolation of optics on rh
         494             :    real(r8), intent(in) :: wrh(pcols,pver)  ! weight for linear interpolation of optics on rh
         495             :    real(r8), intent(in) :: mass(pcols,pver)
         496             :    real(r8), intent(in) :: ext(:,:)
         497             :    real(r8), intent(in) :: ssa(:,:)
         498             :    real(r8), intent(in) :: asm(:,:)
         499             :
         500             :    real(r8), intent(out) :: tau    (pcols,pver,nswbands)
         501             :    real(r8), intent(out) :: tau_w  (pcols,pver,nswbands)
         502             :    real(r8), intent(out) :: tau_w_g(pcols,pver,nswbands)
         503             :    real(r8), intent(out) :: tau_w_f(pcols,pver,nswbands)
         504             :
         505             :    ! Local variables
         506             :    real(r8) :: ext1, ssa1, asm1
         507             :    integer :: icol, ilev, iswband
         508             :    !-----------------------------------------------------------------------------
         509             :
         510    67152240 :    do iswband = 1, nswbands
         511  1051010640 :       do icol = 1, ncol
         512 92545365024 :          do ilev = 1, pver
         513 >36599*10^7 :             ext1 = (1 + wrh(icol,ilev)) * ext(krh(icol,ilev)+1,iswband) &
         514 >45749*10^7 :                       - wrh(icol,ilev)  * ext(krh(icol,ilev),  iswband)
         515 91498831200 :             ssa1 = (1 + wrh(icol,ilev)) * ssa(krh(icol,ilev)+1,iswband) &
         516 91498831200 :                       - wrh(icol,ilev)  * ssa(krh(icol,ilev),  iswband)
         517 91498831200 :             asm1 = (1 + wrh(icol,ilev)) * asm(krh(icol,ilev)+1,iswband) &
         518 91498831200 :                       - wrh(icol,ilev)  * asm(krh(icol,ilev),  iswband)
         519             :
         520 91498831200 :             tau    (icol, ilev, iswband) = mass(icol, ilev) * ext1
         521 91498831200 :             tau_w  (icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1
         522 91498831200 :             tau_w_g(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1
         523 92482689600 :             tau_w_f(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1 * asm1
         524             :          enddo
         525             :       enddo
         526             :    enddo
         527             :
         528      746136 : end subroutine get_hygro_rad_props
         529             :
         530             : !==============================================================================
         531             :
         532     5222952 : subroutine get_nonhygro_rad_props(ncol, mass, ext, ssa, asm, &
         533             :                                   tau, tau_w, tau_w_g, tau_w_f)
         534             :
         535             :    ! Arguments
         536             :    integer,  intent(in) :: ncol
         537             :    real(r8), intent(in) :: mass(pcols, pver)
         538             :    real(r8), intent(in) :: ext(:)
         539             :    real(r8), intent(in) :: ssa(:)
         540             :    real(r8), intent(in) :: asm(:)
         541             :
         542             :    real(r8), intent(out) :: tau    (pcols, pver, nswbands)
         543             :    real(r8), intent(out) :: tau_w  (pcols, pver, nswbands)
         544             :    real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
         545             :    real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
         546             :
         547             :    ! Local variables
         548             :    integer  :: iswband
         549             :    real(r8) :: ext1, ssa1, asm1
         550             :    !-----------------------------------------------------------------------------
         551             :
         552    78344280 :    do iswband = 1, nswbands
         553    73121328 :       ext1 = ext(iswband)
         554    73121328 :       ssa1 = ssa(iswband)
         555    73121328 :       asm1 = asm(iswband)
         556 >11362*10^7 :       tau    (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1
         557 >11362*10^7 :       tau_w  (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1
         558 >11362*10^7 :       tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1
         559 >11362*10^7 :       tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1 * asm1
         560             :    enddo
         561             :
         562     5222952 : end subroutine get_nonhygro_rad_props
         563             :
         564             : !==============================================================================
         565             :
         566           0 : subroutine get_volcanic_radius_rad_props(ncol, mass, pbuf_radius_name, pbuf,  r_ext, r_scat, r_ascat, r_mu, &
         567             :                                   tau, tau_w, tau_w_g, tau_w_f)
         568             :
         569             :
         570             :    use physics_buffer, only : pbuf_get_field, pbuf_get_index
         571             :
         572             :    ! Arguments
         573             :    integer,  intent(in) :: ncol
         574             :    real(r8), intent(in) :: mass(pcols, pver)
         575             :    character(len=*) :: pbuf_radius_name
         576             :    type(physics_buffer_desc), pointer :: pbuf(:)
         577             :    real(r8), intent(in) :: r_ext(:,:)
         578             :    real(r8), intent(in) :: r_scat(:,:)
         579             :    real(r8), intent(in) :: r_ascat(:,:)
         580             :    real(r8), intent(in) :: r_mu(:) ! log(radius) domain of mass-specific optics
         581             :
         582             :    real(r8), intent(out) :: tau    (pcols, pver, nswbands)
         583             :    real(r8), intent(out) :: tau_w  (pcols, pver, nswbands)
         584             :    real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
         585             :    real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
         586             :
         587             :    ! Local variables
         588             :    integer  :: iswband
         589             :    real(r8) :: g
         590             :
         591             :    integer  :: idx                             ! index to radius in physics buffer
         592           0 :    real(r8), pointer :: geometric_radius(:,:)  ! geometric mean radius of volcanic aerosol
         593             :    real(r8) :: mu(pcols,pver)                  ! log(geometric mean radius of volcanic aerosol)
         594             :    integer  :: kmu, nmu
         595             :    real(r8) :: wmu, mutrunc, r_mu_max, r_mu_min
         596             :
         597             :    ! interpolated values from table
         598             :    real(r8) :: ext(nswbands)
         599             :    real(r8) :: scat(nswbands)
         600             :    real(r8) :: ascat(nswbands)
         601             :
         602             :    integer :: i, k ! column level iterator
         603             :    !-----------------------------------------------------------------------------
         604             :
         605           0 :    tau    =0._r8
         606           0 :    tau_w  =0._r8
         607           0 :    tau_w_g=0._r8
         608           0 :    tau_w_f=0._r8
         609             :
         610             :    ! get microphysical properties for volcanic aerosols
         611           0 :    idx = pbuf_get_index(pbuf_radius_name)
         612           0 :    call pbuf_get_field(pbuf, idx, geometric_radius )
         613             :
         614             :    ! interpolate in radius
         615             :    ! caution: clip the table with no warning when outside bounds
         616           0 :    nmu = size(r_mu)
         617           0 :    r_mu_max = r_mu(nmu)
         618           0 :    r_mu_min = r_mu(1)
         619           0 :    do i = 1, ncol
         620           0 :       do k = 1, pver
         621           0 :          if(geometric_radius(i,k) > 0._r8) then
         622           0 :             mu(i,k) = log(geometric_radius(i,k))
         623             :          else
         624           0 :             mu(i,k) = 0._r8
         625             :          endif
         626           0 :          mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min)
         627           0 :          kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8)
         628           0 :          wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8)
         629           0 :          do iswband = 1, nswbands
         630           0 :             ext(iswband) =  &
         631           0 :                ((1._r8 - wmu) * r_ext(iswband, kmu  ) + &
         632           0 :                (wmu) * r_ext(iswband, kmu+1))
         633             :             scat(iswband) =  &
         634           0 :                ((1._r8 - wmu) * r_scat(iswband, kmu  ) + &
         635           0 :                (wmu) * r_scat(iswband, kmu+1))
         636             :             ascat(iswband) =  &
         637           0 :                ((1._r8 - wmu) * r_ascat(iswband, kmu  ) + &
         638           0 :                (wmu) * r_ascat(iswband, kmu+1))
         639           0 :             if (scat(iswband).gt.0._r8) then
         640           0 :                g = ascat(iswband)/scat(iswband)
         641             :             else
         642             :                g=0._r8
         643             :             endif
         644           0 :             tau    (i,k,iswband) = mass(i,k) * ext(iswband)
         645           0 :             tau_w  (i,k,iswband) = mass(i,k) * scat(iswband)
         646           0 :             tau_w_g(i,k,iswband) = mass(i,k) * ascat(iswband)
         647           0 :             tau_w_f(i,k,iswband) = mass(i,k) * g * ascat(iswband)
         648             :          end do
         649             :       enddo
         650             :    enddo
         651             :
         652           0 : end subroutine get_volcanic_radius_rad_props
         653             :
         654             : !==============================================================================
         655             :
         656           0 : subroutine get_volcanic_rad_props(ncol, mass, ext, scat, ascat, &
         657             :                                   tau, tau_w, tau_w_g, tau_w_f)
         658             :
         659             :    ! Arguments
         660             :    integer,  intent(in) :: ncol
         661             :    real(r8), intent(in) :: mass(pcols, pver)
         662             :    real(r8), intent(in) :: ext(:)
         663             :    real(r8), intent(in) :: scat(:)
         664             :    real(r8), intent(in) :: ascat(:)
         665             :
         666             :    real(r8), intent(out) :: tau    (pcols, pver, nswbands)
         667             :    real(r8), intent(out) :: tau_w  (pcols, pver, nswbands)
         668             :    real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
         669             :    real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
         670             :
         671             :    ! Local variables
         672             :    integer  :: iswband
         673             :    real(r8) :: g
         674             :    !-----------------------------------------------------------------------------
         675             :
         676           0 :    do iswband = 1, nswbands
         677           0 :       if (scat(iswband).gt.0._r8) then
         678           0 :          g = ascat(iswband)/scat(iswband)
         679             :       else
         680             :          g=0._r8
         681             :       endif
         682           0 :       tau    (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext(iswband)
         683           0 :       tau_w  (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * scat(iswband)
         684           0 :       tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ascat(iswband)
         685           0 :       tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * g * ascat(iswband)
         686             :    enddo
         687             :
         688           0 : end subroutine get_volcanic_rad_props
         689             :
         690             : !==============================================================================
         691             :
         692    10445904 : subroutine aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaer, tau, diag_idx, troplev)
         693             :
         694             :    ! output aerosol optical depth for the visible band
         695             :
         696             :    integer,          intent(in) :: lchnk
         697             :    integer,          intent(in) :: ncol           ! number of columns
         698             :    integer,          intent(in) :: nnite          ! number of night columns
         699             :    integer,          intent(in) :: idxnite(:)     ! local column indices of night columns
         700             :    integer,          intent(in) :: iaer           ! aerosol index -- if 0 then tau is a total for all aerosols
         701             :    real(r8),         intent(in) :: tau(:,:)       ! aerosol optical depth for the visible band
         702             :    integer,          intent(in) :: diag_idx       ! identifies whether the aerosol optics
         703             :                                                   ! is for the climate calc or a diagnostic calc
         704             :    integer,          intent(in) :: troplev(:)     ! tropopause level
         705             :
         706             :    ! Local variables
         707             :    integer  :: i
         708             :    real(r8) :: tmp(pcols), tmp2(pcols)
         709             :    !-----------------------------------------------------------------------------
         710             :
         711             :    ! currently only implemented for climate calc
         712    10445904 :    if (diag_idx > 0) return
         713             :
         714             :    ! compute total column aerosol optical depth
         715 15435940104 :    tmp(1:ncol) = sum(tau(1:ncol,:), 2)
         716             :    ! use fillvalue to indicate night columns
         717    92434104 :    do i = 1, nnite
         718    92434104 :       tmp(idxnite(i)) = fillvalue
         719             :    end do
         720             :
         721    10445904 :    if (iaer > 0) then
         722     9699768 :       call outfld(odv_names(iaer), tmp, pcols, lchnk)
         723             :    else
         724      746136 :       call outfld('AEROD_v', tmp, pcols, lchnk)
         725    12458736 :       do i = 1, ncol
         726   624127728 :         tmp2(i) = sum(tau(i,:troplev(i)))
         727             :       end do
         728      746136 :       call outfld('AODvstrt', tmp2, pcols, lchnk)
         729             :    end if
         730             :
         731             : end subroutine aer_vis_diag_out
         732             :
         733             : !==============================================================================
         734             :
         735             : end module aer_rad_props

Generated by: LCOV version 1.14