LCOV - code coverage report
Current view: | top level - physics/cam - aerosol_optics_cam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 11 | 667 | 1.6 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 14 | 7.1 %

              Line data    Source code
           1             : module aerosol_optics_cam
           2             :   use shr_kind_mod, only: r8 => shr_kind_r8
           3             :   use shr_kind_mod, only: cl => shr_kind_cl
           4             :   use cam_logfile,  only: iulog
           5             :   use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_uv_diag, idx_nir_diag
           6             :   use radconstants, only: get_lw_spectral_boundaries
           7             :   use phys_prop,    only: ot_length
           8             :   use physics_types,only: physics_state
           9             :   use physics_buffer,only: physics_buffer_desc
          10             :   use ppgrid, only: pcols, pver
          11             :   use physconst, only: rga, rair
          12             :   use cam_abortutils, only: endrun
          13             :   use spmd_utils, only: masterproc
          14             :   use rad_constituents,  only: n_diag, rad_cnst_get_call_list
          15             :   use cam_history,       only: addfld, add_default, outfld, horiz_only, fieldname_len
          16             :   use cam_history_support, only: fillvalue
          17             :
          18             :   use tropopause, only : tropopause_findChemTrop
          19             :
          20             :   use aerosol_properties_mod, only: aerosol_properties
          21             :   use modal_aerosol_properties_mod, only: modal_aerosol_properties
          22             :
          23             :   use aerosol_state_mod,      only: aerosol_state
          24             :   use modal_aerosol_state_mod,only: modal_aerosol_state
          25             :
          26             :   use aerosol_optics_mod,     only: aerosol_optics
          27             :   use refractive_aerosol_optics_mod, only: refractive_aerosol_optics
          28             :
          29             :   implicit none
          30             :
          31             :   private
          32             :
          33             :   public :: aerosol_optics_cam_readnl
          34             :   public :: aerosol_optics_cam_init
          35             :   public :: aerosol_optics_cam_final
          36             :   public :: aerosol_optics_cam_sw
          37             :   public :: aerosol_optics_cam_lw
          38             :
          39             :   type aero_props_t
          40             :      class(aerosol_properties), pointer :: obj => null()
          41             :   end type aero_props_t
          42             :   type aero_state_t
          43             :      class(aerosol_state), pointer :: obj => null()
          44             :   end type aero_state_t
          45             :
          46             :   type(aero_props_t), allocatable :: aero_props(:) ! array of aerosol properties objects to allow for
          47             :                                                    ! multiple aerosol representations in the same sim
          48             :                                                    ! such as MAM and CARMA
          49             :
          50             :   ! refractive index for water read in read_water_refindex
          51             :   complex(r8) :: crefwsw(nswbands) = -huge(1._r8) ! complex refractive index for water visible
          52             :   complex(r8) :: crefwlw(nlwbands) = -huge(1._r8) ! complex refractive index for water infrared
          53             :   character(len=cl) :: water_refindex_file = 'NONE' ! full pathname for water refractive index dataset
          54             :
          55             :   logical :: modal_active = .false.
          56             :   integer :: num_aero_models = 0
          57             :   integer :: lw10um_indx = -1            ! wavelength index corresponding to 10 microns
          58             :   real(r8), parameter :: lw10um = 10._r8 ! microns
          59             :
          60             :   character(len=4) :: diag(0:n_diag) = (/'    ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ', '_d6 ','_d7 ','_d8 ','_d9 ','_d10'/)
          61             :
          62             :   type out_name
          63             :      character(len=fieldname_len), allocatable :: name(:) ! nbins
          64             :   end type out_name
          65             :
          66             :   type(out_name), allocatable :: burden_fields(:) ! num_aero_models
          67             :   type(out_name), allocatable :: aodbin_fields(:)
          68             :   type(out_name), allocatable :: aoddust_fields(:)
          69             :   type(out_name), allocatable :: burdendn_fields(:) ! num_aero_models
          70             :   type(out_name), allocatable :: aodbindn_fields(:)
          71             :   type(out_name), allocatable :: aoddustdn_fields(:)
          72             :
          73             : contains
          74             :
          75             :   !===============================================================================
          76        1536 :   subroutine aerosol_optics_cam_readnl(nlfile)
          77             :     use namelist_utils, only : find_group_name
          78             :     use spmd_utils,     only : mpicom, masterprocid, mpi_character, mpi_success
          79             :
          80             :     character(len=*), intent(in)  :: nlfile  ! filepath for file containing namelist input
          81             :
          82             :     integer                       :: unitn, ierr
          83             :     character(len=cl)             :: errmsg
          84             :     character(len=*), parameter   :: subname = 'aerosol_optics_cam_readnl'
          85             :
          86             :     ! ===================
          87             :     ! Namelist definition
          88             :     ! ===================
          89             :     namelist /aerosol_optics_nl/ water_refindex_file
          90             :
          91             :     ! =============
          92             :     ! Read namelist
          93             :     ! =============
          94        1536 :     if (masterproc) then
          95           2 :        open( newunit=unitn, file=trim(nlfile), status='old' )
          96           2 :        call find_group_name(unitn, 'aerosol_optics_nl', status=ierr)
          97           2 :        if (ierr == 0) then
          98           0 :           read(unitn, aerosol_optics_nl, iostat=ierr)
          99           0 :           if (ierr /= 0) then
         100           0 :              write(errmsg,'(2a,i10)') subname,':: ERROR reading namelist, error code: ',ierr
         101           0 :              call endrun(errmsg)
         102             :           end if
         103             :        end if
         104           2 :        close(unitn)
         105             :     end if
         106             :
         107             :     ! ============================
         108             :     ! Broadcast namelist variables
         109             :     ! ============================
         110        1536 :     call mpi_bcast(water_refindex_file, len(water_refindex_file),  mpi_character, masterprocid, mpicom, ierr)
         111        1536 :     if (ierr/=mpi_success) then
         112           0 :        call endrun(subname // ':: ERROR mpi_bcast '//trim(water_refindex_file))
         113             :     end if
         114             :
         115        1536 :     if (masterproc) then
         116           2 :        write(iulog,*) subname,': water_refindex_file = ',trim(water_refindex_file)
         117             :     end if
         118             :
         119        1536 :   end subroutine aerosol_optics_cam_readnl
         120             :
         121             :   !===============================================================================
         122           0 :   subroutine aerosol_optics_cam_init
         123             :     use rad_constituents, only: rad_cnst_get_info
         124             :     use phys_control,     only: phys_getopts
         125             :     use ioFileMod,        only: getfil
         126             :
         127             :     character(len=*), parameter :: prefix = 'aerosol_optics_cam_init: '
         128             :     integer :: nmodes=0, iaermod, istat, ilist, i
         129             :
         130             :     logical :: call_list(0:n_diag)
         131             :     real(r8) :: lwavlen_lo(nlwbands), lwavlen_hi(nlwbands)
         132             :     integer :: m, n
         133             :
         134             :     character(len=fieldname_len) :: fldname
         135             :     character(len=128) :: lngname
         136             :     logical :: history_aero_optics     ! output aerosol optics diagnostics
         137             :     logical :: history_amwg            ! output the variables used by the AMWG diag package
         138             :     logical :: history_dust            ! output dust diagnostics
         139             :
         140             :     character(len=cl) :: locfile
         141             :
         142             :     call phys_getopts(history_amwg_out        = history_amwg, &
         143             :                       history_aero_optics_out = history_aero_optics, &
         144           0 :                       history_dust_out        = history_dust )
         145             :
         146           0 :     num_aero_models = 0
         147             :
         148           0 :     call rad_cnst_get_info(0, nmodes=nmodes)
         149           0 :     modal_active = nmodes>0
         150             :
         151           0 :     if (modal_active) then
         152           0 :        num_aero_models = num_aero_models+1 ! count aerosol models
         153             :     end if
         154             :
         155           0 :     if (num_aero_models>0) then
         156           0 :        allocate(aero_props(num_aero_models), stat=istat)
         157           0 :        if (istat/=0) then
         158           0 :           call endrun(prefix//'array allocation error: aero_props')
         159             :        end if
         160             :     end if
         161             :
         162           0 :     iaermod = 0
         163             :
         164           0 :     if (modal_active) then
         165           0 :        iaermod = iaermod+1
         166           0 :        aero_props(iaermod)%obj => modal_aerosol_properties()
         167             :     end if
         168             :
         169           0 :     if (water_refindex_file=='NONE') then
         170           0 :        call endrun(prefix//'water_refindex_file must be specified')
         171             :     else
         172           0 :        call getfil(water_refindex_file, locfile)
         173           0 :        call read_water_refindex(locfile)
         174             :     end if
         175             :
         176           0 :     call get_lw_spectral_boundaries(lwavlen_lo, lwavlen_hi, units='um')
         177           0 :     do i = 1,nlwbands
         178           0 :        if ((lwavlen_lo(i)<=lw10um) .and. (lwavlen_hi(i)>=lw10um)) then
         179           0 :           lw10um_indx = i ! index corresponding to 10 microns
         180             :        end if
         181             :     end do
         182           0 :     call rad_cnst_get_call_list(call_list)
         183             :
         184           0 :     do ilist = 0, n_diag
         185           0 :        if (call_list(ilist)) then
         186             :           call addfld ('EXTINCT'//diag(ilist),    (/ 'lev' /), 'A','/m',&
         187           0 :                'Aerosol extinction 550 nm, day only', flag_xyfill=.true.)
         188             :           call addfld ('EXTINCTUV'//diag(ilist),  (/ 'lev' /), 'A','/m',&
         189           0 :                'Aerosol extinction 350 nm, day only', flag_xyfill=.true.)
         190             :           call addfld ('EXTINCTNIR'//diag(ilist), (/ 'lev' /), 'A','/m',&
         191           0 :                'Aerosol extinction 1020 nm, day only', flag_xyfill=.true.)
         192             :           call addfld ('ABSORB'//diag(ilist),     (/ 'lev' /), 'A','/m',&
         193           0 :                'Aerosol absorption, day only', flag_xyfill=.true.)
         194             :           call addfld ('AODVIS'//diag(ilist),   horiz_only,  'A','  ', &
         195           0 :                'Aerosol optical depth 550 nm', flag_xyfill=.true.)
         196             :           call addfld ('AODVISst'//diag(ilist), horiz_only,  'A','  ', &
         197           0 :                'Stratospheric aerosol optical depth 550 nm, day only', flag_xyfill=.true.)
         198             :           call addfld ('AODNIRst'//diag(ilist), horiz_only,  'A','  ', &
         199           0 :                'Stratospheric aerosol optical depth 1020 nm, day only',flag_xyfill=.true.)
         200             :           call addfld ('AODUVst'//diag(ilist),  horiz_only,  'A','  ', &
         201           0 :                'Stratospheric aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
         202             :           call addfld ('AODUV'//diag(ilist),      horiz_only,  'A','  ', &
         203           0 :                'Aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
         204             :           call addfld ('AODNIR'//diag(ilist),     horiz_only,  'A','  ', &
         205           0 :                'Aerosol optical depth 1020 nm, day only',flag_xyfill=.true.)
         206             :           call addfld ('AODABS'//diag(ilist),     horiz_only,  'A','  ', &
         207           0 :                'Aerosol absorption optical depth 550 nm, day only', flag_xyfill=.true.)
         208             :           call addfld ('AODxASYM'//diag(ilist),   horiz_only,  'A','  ', &
         209           0 :                'Aerosol optical depth 550 * asymmetry factor, day only', flag_xyfill=.true.)
         210             :           call addfld ('EXTxASYM'//diag(ilist),   (/ 'lev' /), 'A','  ', &
         211           0 :                'extinction 550 nm * asymmetry factor, day only',  flag_xyfill=.true.)
         212             :           call addfld ('AODTOT'//diag(ilist), horiz_only, 'A','1',&
         213           0 :                'Aerosol optical depth summed over all sw wavelengths', flag_xyfill=.true.)
         214             :
         215             :           call addfld ('EXTINCTdn'//diag(ilist),    (/ 'lev' /), 'A','/m',&
         216           0 :                'Aerosol extinction 550 nm, day only', flag_xyfill=.true.)
         217             :           call addfld ('EXTINCTUVdn'//diag(ilist),  (/ 'lev' /), 'A','/m',&
         218           0 :                'Aerosol extinction 350 nm, day only', flag_xyfill=.true.)
         219             :           call addfld ('EXTINCTNIRdn'//diag(ilist), (/ 'lev' /), 'A','/m',&
         220           0 :                'Aerosol extinction 1020 nm, day only', flag_xyfill=.true.)
         221             :           call addfld ('ABSORBdn'//diag(ilist),     (/ 'lev' /), 'A','/m',&
         222           0 :                'Aerosol absorption, day only', flag_xyfill=.true.)
         223             :           call addfld ('AODVISdn'//diag(ilist),   horiz_only,  'A','  ', &
         224           0 :                'Aerosol optical depth 550 nm', flag_xyfill=.true.)
         225             :           call addfld ('AODVISstdn'//diag(ilist), horiz_only,  'A','  ', &
         226           0 :                'Stratospheric aerosol optical depth 550 nm, day only', flag_xyfill=.true.)
         227             :           call addfld ('AODNIRstdn'//diag(ilist), horiz_only,  'A','  ', &
         228           0 :                'Stratospheric aerosol optical depth 1020 nm, day only', flag_xyfill=.true.)
         229             :           call addfld ('AODUVstdn'//diag(ilist),  horiz_only,  'A','  ', &
         230           0 :                'Stratospheric aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
         231             :           call addfld ('AODUVdn'//diag(ilist),      horiz_only,  'A','  ', &
         232           0 :                'Aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
         233             :           call addfld ('AODNIRdn'//diag(ilist),     horiz_only,  'A','  ', &
         234           0 :                'Aerosol optical depth 1020 nm, day only', flag_xyfill=.true.)
         235             :           call addfld ('AODABSdn'//diag(ilist),     horiz_only,  'A','  ', &
         236           0 :                'Aerosol absorption optical depth 550 nm, day only', flag_xyfill=.true.)
         237             :           call addfld ('AODxASYMdn'//diag(ilist),   horiz_only,  'A','  ', &
         238           0 :                'Aerosol optical depth 550 * asymmetry factor, day only', flag_xyfill=.true.)
         239             :           call addfld ('EXTxASYMdn'//diag(ilist),   (/ 'lev' /), 'A','  ', &
         240           0 :                'extinction 550 nm * asymmetry factor, day only',  flag_xyfill=.true.)
         241             :           call addfld ('AODTOTdn'//diag(ilist), horiz_only, 'A','1',&
         242           0 :                'Aerosol optical depth summed over all sw wavelengths, day only')
         243             :
         244           0 :           if (lw10um_indx>0) then
         245             :              call addfld('AODABSLW'//diag(ilist), (/ 'lev' /), 'A','/m',&
         246           0 :                   'Aerosol long-wave absorption optical depth at 10 microns')
         247             :           end if
         248             :           call addfld ('TOTABSLW'//diag(ilist), (/ 'lev' /), 'A',' ', &
         249           0 :                'LW Aero total abs')
         250             :
         251           0 :           if (ilist>0 .and. history_aero_optics) then
         252           0 :              call add_default ('EXTINCT'//diag(ilist), 1, ' ')
         253           0 :              call add_default ('ABSORB'//diag(ilist),  1, ' ')
         254           0 :              call add_default ('AODVIS'//diag(ilist),  1, ' ')
         255           0 :              call add_default ('AODVISst'//diag(ilist),  1, ' ')
         256           0 :              call add_default ('AODABS'//diag(ilist),  1, ' ')
         257             :           end if
         258             :
         259             :        end if
         260             :     end do
         261             :
         262           0 :     if (num_aero_models>0) then
         263             :
         264           0 :        allocate(burden_fields(num_aero_models), stat=istat)
         265           0 :        if (istat/=0) then
         266           0 :           call endrun(prefix//'array allocation error: burden_fields')
         267             :        end if
         268           0 :        allocate(aodbin_fields(num_aero_models), stat=istat)
         269           0 :        if (istat/=0) then
         270           0 :           call endrun(prefix//'array allocation error: aodbin_fields')
         271             :        end if
         272           0 :        allocate(aoddust_fields(num_aero_models), stat=istat)
         273           0 :        if (istat/=0) then
         274           0 :           call endrun(prefix//'array allocation error: aoddust_fields')
         275             :        end if
         276             :
         277           0 :        allocate(burdendn_fields(num_aero_models), stat=istat)
         278           0 :        if (istat/=0) then
         279           0 :           call endrun(prefix//'array allocation error: burdendn_fields')
         280             :        end if
         281           0 :        allocate(aodbindn_fields(num_aero_models), stat=istat)
         282           0 :        if (istat/=0) then
         283           0 :           call endrun(prefix//'array allocation error: aodbindn_fields')
         284             :        end if
         285           0 :        allocate(aoddustdn_fields(num_aero_models), stat=istat)
         286           0 :        if (istat/=0) then
         287           0 :           call endrun(prefix//'array allocation error: aoddustdn_fields')
         288             :        end if
         289             :
         290           0 :        do n = 1,num_aero_models
         291             :
         292           0 :           allocate(burden_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         293           0 :           if (istat/=0) then
         294           0 :              call endrun(prefix//'array allocation error: burden_fields(n)%name')
         295             :           end if
         296           0 :           allocate(aodbin_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         297           0 :           if (istat/=0) then
         298           0 :              call endrun(prefix//'array allocation error: aodbin_fields(n)%name')
         299             :           end if
         300           0 :           allocate(aoddust_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         301           0 :           if (istat/=0) then
         302           0 :              call endrun(prefix//'array allocation error: aoddust_fields(n)%name')
         303             :           end if
         304             :
         305           0 :           allocate(burdendn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         306           0 :           if (istat/=0) then
         307           0 :              call endrun(prefix//'array allocation error: burdendn_fields(n)%name')
         308             :           end if
         309           0 :           allocate(aodbindn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         310           0 :           if (istat/=0) then
         311           0 :              call endrun(prefix//'array allocation error: aodbindn_fields(n)%name')
         312             :           end if
         313           0 :           allocate(aoddustdn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
         314           0 :           if (istat/=0) then
         315           0 :              call endrun(prefix//'array allocation error: aoddustdn_fields(n)%name')
         316             :           end if
         317             :
         318           0 :           do m = 1, aero_props(n)%obj%nbins()
         319             :
         320           0 :              write(fldname,'(a,i2.2)') 'BURDEN', m
         321           0 :              burden_fields(n)%name(m) = fldname
         322           0 :              write(lngname,'(a,i2.2)') 'Aerosol burden bin ', m
         323           0 :              call addfld (fldname, horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.)
         324           0 :              if (history_aero_optics) then
         325           0 :                 call add_default (fldname, 1, ' ')
         326             :              end if
         327             :
         328           0 :              fldname = 'AOD_'//trim(aero_props(n)%obj%bin_name(0,m))
         329           0 :              aodbin_fields(n)%name(m) = fldname
         330           0 :              lngname = 'Aerosol optical depth, day only, 550 nm, '//trim(aero_props(n)%obj%bin_name(0,m))
         331           0 :              call addfld (aodbin_fields(n)%name(m), horiz_only, 'A', '  ', lngname, flag_xyfill=.true.)
         332           0 :              if (history_aero_optics) then
         333           0 :                 call add_default (fldname, 1, ' ')
         334             :              end if
         335             :
         336           0 :              write(fldname,'(a,i2.2)') 'AODDUST', m
         337           0 :              aoddust_fields(n)%name(m) = fldname
         338           0 :              write(lngname,'(a,i2,a)') 'Aerosol optical depth, day only, 550 nm mode ',m,' from dust'
         339           0 :              call addfld (aoddust_fields(n)%name(m), horiz_only, 'A', '  ', lngname, flag_xyfill=.true.)
         340           0 :              if (history_aero_optics) then
         341           0 :                 call add_default (fldname, 1, ' ')
         342             :              end if
         343             :
         344           0 :              write(fldname,'(a,i2.2)') 'BURDENdn', m
         345           0 :              burdendn_fields(n)%name(m) = fldname
         346           0 :              write(lngname,'(a,i2)') 'Aerosol burden, day night, bin ', m
         347           0 :              call addfld (burdendn_fields(n)%name(m), horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.)
         348           0 :              if (history_aero_optics) then
         349           0 :                 call add_default (fldname, 1, ' ')
         350             :              end if
         351             :
         352           0 :              fldname = 'AODdn_'//trim(aero_props(n)%obj%bin_name(0,m))
         353           0 :              aodbindn_fields(n)%name(m) = fldname
         354           0 :              lngname = 'Aerosol optical depth 550 nm, day night, '//trim(aero_props(n)%obj%bin_name(0,m))
         355           0 :              call addfld (aodbindn_fields(n)%name(m), horiz_only, 'A', '  ', lngname, flag_xyfill=.true.)
         356           0 :              if (history_aero_optics) then
         357           0 :                 call add_default (fldname, 1, ' ')
         358             :              end if
         359             :
         360           0 :              write(fldname,'(a,i2.2)') 'AODdnDUST', m
         361           0 :              aoddustdn_fields(n)%name(m) = fldname
         362           0 :              write(lngname,'(a,i2,a)') 'Aerosol optical depth 550 nm, day night, bin ',m,' from dust'
         363           0 :              call addfld (aoddustdn_fields(n)%name(m), horiz_only, 'A', '  ', lngname, flag_xyfill=.true.)
         364           0 :              if (history_aero_optics) then
         365           0 :                 call add_default (fldname, 1, ' ')
         366             :              end if
         367             :
         368             :           end do
         369             :
         370             :        end do
         371             :
         372             :     end if
         373             :
         374             :    call addfld ('AODDUST',       horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from dust, day only',         &
         375           0 :         flag_xyfill=.true.)
         376             :    call addfld ('AODSO4',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from SO4, day only',          &
         377           0 :         flag_xyfill=.true.)
         378             :    call addfld ('AODPOM',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from POM, day only',          &
         379           0 :         flag_xyfill=.true.)
         380             :    call addfld ('AODSOA',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from SOA, day only',          &
         381           0 :         flag_xyfill=.true.)
         382             :    call addfld ('AODBC',         horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from BC, day only',           &
         383           0 :         flag_xyfill=.true.)
         384             :    call addfld ('AODSS',         horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from seasalt, day only',      &
         385           0 :         flag_xyfill=.true.)
         386             :    call addfld ('AODABSBC',      horiz_only, 'A','  ',    'Aerosol absorption optical depth 550 nm from BC, day only',&
         387           0 :         flag_xyfill=.true.)
         388             :    call addfld ('BURDENDUST',    horiz_only, 'A','kg/m2', 'Dust aerosol burden, day only'        ,                    &
         389           0 :         flag_xyfill=.true.)
         390             :    call addfld ('BURDENSO4',     horiz_only, 'A','kg/m2', 'Sulfate aerosol burden, day only'     ,                    &
         391           0 :         flag_xyfill=.true.)
         392             :    call addfld ('BURDENPOM',     horiz_only, 'A','kg/m2', 'POM aerosol burden, day only'         ,                    &
         393           0 :         flag_xyfill=.true.)
         394             :    call addfld ('BURDENSOA',     horiz_only, 'A','kg/m2', 'SOA aerosol burden, day only'         ,                    &
         395           0 :         flag_xyfill=.true.)
         396             :    call addfld ('BURDENBC',      horiz_only, 'A','kg/m2', 'Black carbon aerosol burden, day only',                    &
         397           0 :         flag_xyfill=.true.)
         398             :    call addfld ('BURDENSEASALT', horiz_only, 'A','kg/m2', 'Seasalt aerosol burden, day only'     ,                    &
         399           0 :         flag_xyfill=.true.)
         400             :    call addfld ('SSAVIS',        horiz_only, 'A','  ',    'Aerosol single-scatter albedo, day only',                  &
         401           0 :         flag_xyfill=.true.)
         402             :
         403             :    call addfld ('AODDUSTdn',       horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from dust, day night',         &
         404           0 :         flag_xyfill=.true.)
         405             :    call addfld ('AODSO4dn',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from SO4, day night',          &
         406           0 :         flag_xyfill=.true.)
         407             :    call addfld ('AODPOMdn',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from POM, day night',          &
         408           0 :         flag_xyfill=.true.)
         409             :    call addfld ('AODSOAdn',        horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from SOA, day night',          &
         410           0 :         flag_xyfill=.true.)
         411             :    call addfld ('AODBCdn',         horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from BC, day night',           &
         412           0 :         flag_xyfill=.true.)
         413             :    call addfld ('AODSSdn',         horiz_only, 'A','  ',    'Aerosol optical depth 550 nm from seasalt, day night',      &
         414           0 :         flag_xyfill=.true.)
         415             :    call addfld ('AODABSBCdn',      horiz_only, 'A','  ',    'Aerosol absorption optical depth 550 nm from BC, day night',&
         416           0 :         flag_xyfill=.true.)
         417             :    call addfld ('BURDENDUSTdn',    horiz_only, 'A','kg/m2', 'Dust aerosol burden, day night'        ,                    &
         418           0 :         flag_xyfill=.true.)
         419             :    call addfld ('BURDENSO4dn',     horiz_only, 'A','kg/m2', 'Sulfate aerosol burden, day night'     ,                    &
         420           0 :         flag_xyfill=.true.)
         421             :    call addfld ('BURDENPOMdn',     horiz_only, 'A','kg/m2', 'POM aerosol burden, day night'         ,                    &
         422           0 :         flag_xyfill=.true.)
         423             :    call addfld ('BURDENSOAdn',     horiz_only, 'A','kg/m2', 'SOA aerosol burden, day night'         ,                    &
         424           0 :         flag_xyfill=.true.)
         425             :    call addfld ('BURDENBCdn',      horiz_only, 'A','kg/m2', 'Black carbon aerosol burden, day night',                    &
         426           0 :         flag_xyfill=.true.)
         427             :    call addfld ('BURDENSEASALTdn', horiz_only, 'A','kg/m2', 'Seasalt aerosol burden, day night'     ,                    &
         428           0 :         flag_xyfill=.true.)
         429             :    call addfld ('SSAVISdn',        horiz_only, 'A','  ',    'Aerosol single-scatter albedo, day night',                  &
         430           0 :         flag_xyfill=.true.)
         431             :
         432           0 :    if (history_amwg) then
         433           0 :       call add_default ('AODDUST01'     , 1, ' ')
         434           0 :       call add_default ('AODDUST03'     , 1, ' ')
         435           0 :       call add_default ('AODDUST'      , 1, ' ')
         436           0 :       call add_default ('AODVIS'       , 1, ' ')
         437             :    end if
         438             :
         439           0 :    if (history_dust) then
         440           0 :       call add_default ('AODDUST01'     , 1, ' ')
         441           0 :       call add_default ('AODDUST02'     , 1, ' ')
         442           0 :       call add_default ('AODDUST03'     , 1, ' ')
         443             :    end if
         444             :
         445           0 :    if (history_aero_optics) then
         446           0 :       call add_default ('AODDUST01'     , 1, ' ')
         447           0 :       call add_default ('AODDUST03'     , 1, ' ')
         448           0 :       call add_default ('ABSORB'       , 1, ' ')
         449           0 :       call add_default ('AODVIS'       , 1, ' ')
         450           0 :       call add_default ('AODUV'        , 1, ' ')
         451           0 :       call add_default ('AODNIR'       , 1, ' ')
         452           0 :       call add_default ('AODABS'       , 1, ' ')
         453           0 :       call add_default ('AODABSBC'     , 1, ' ')
         454           0 :       call add_default ('AODDUST'      , 1, ' ')
         455           0 :       call add_default ('AODSO4'       , 1, ' ')
         456           0 :       call add_default ('AODPOM'       , 1, ' ')
         457           0 :       call add_default ('AODSOA'       , 1, ' ')
         458           0 :       call add_default ('AODBC'        , 1, ' ')
         459           0 :       call add_default ('AODSS'        , 1, ' ')
         460           0 :       call add_default ('BURDEN01'      , 1, ' ')
         461           0 :       call add_default ('BURDEN02'      , 1, ' ')
         462           0 :       call add_default ('BURDEN03'      , 1, ' ')
         463           0 :       call add_default ('BURDENDUST'   , 1, ' ')
         464           0 :       call add_default ('BURDENSO4'    , 1, ' ')
         465           0 :       call add_default ('BURDENPOM'    , 1, ' ')
         466           0 :       call add_default ('BURDENSOA'    , 1, ' ')
         467           0 :       call add_default ('BURDENBC'     , 1, ' ')
         468           0 :       call add_default ('BURDENSEASALT', 1, ' ')
         469           0 :       call add_default ('SSAVIS'       , 1, ' ')
         470           0 :       call add_default ('EXTINCT'      , 1, ' ')
         471           0 :       call add_default ('AODxASYM'     , 1, ' ')
         472           0 :       call add_default ('EXTxASYM'     , 1, ' ')
         473             :
         474           0 :       call add_default ('AODdnDUST01'     , 1, ' ')
         475           0 :       call add_default ('AODdnDUST03'     , 1, ' ')
         476           0 :       call add_default ('ABSORBdn'       , 1, ' ')
         477           0 :       call add_default ('AODVISdn'       , 1, ' ')
         478           0 :       call add_default ('AODUVdn'        , 1, ' ')
         479           0 :       call add_default ('AODNIRdn'       , 1, ' ')
         480           0 :       call add_default ('AODABSdn'       , 1, ' ')
         481           0 :       call add_default ('AODABSBCdn'     , 1, ' ')
         482           0 :       call add_default ('AODDUSTdn'      , 1, ' ')
         483           0 :       call add_default ('AODSO4dn'       , 1, ' ')
         484           0 :       call add_default ('AODPOMdn'       , 1, ' ')
         485           0 :       call add_default ('AODSOAdn'       , 1, ' ')
         486           0 :       call add_default ('AODBCdn'        , 1, ' ')
         487           0 :       call add_default ('AODSSdn'        , 1, ' ')
         488           0 :       call add_default ('BURDENdn01'      , 1, ' ')
         489           0 :       call add_default ('BURDENdn02'      , 1, ' ')
         490           0 :       call add_default ('BURDENdn03'      , 1, ' ')
         491           0 :       call add_default ('BURDENDUSTdn'   , 1, ' ')
         492           0 :       call add_default ('BURDENSO4dn'    , 1, ' ')
         493           0 :       call add_default ('BURDENPOMdn'    , 1, ' ')
         494           0 :       call add_default ('BURDENSOAdn'    , 1, ' ')
         495           0 :       call add_default ('BURDENBCdn'     , 1, ' ')
         496           0 :       call add_default ('BURDENSEASALTdn', 1, ' ')
         497           0 :       call add_default ('SSAVISdn'       , 1, ' ')
         498           0 :       call add_default ('EXTINCTdn'      , 1, ' ')
         499           0 :       call add_default ('AODxASYMdn'     , 1, ' ')
         500           0 :       call add_default ('EXTxASYMdn'     , 1, ' ')
         501             :    end if
         502             :
         503           0 :   end subroutine aerosol_optics_cam_init
         504             :
         505             :   !===============================================================================
         506           0 :   subroutine aerosol_optics_cam_final
         507             :
         508             :     integer :: iaermod
         509             :
         510           0 :     do iaermod = 1,num_aero_models
         511           0 :        if (associated(aero_props(iaermod)%obj)) then
         512           0 :           deallocate(aero_props(iaermod)%obj)
         513           0 :           nullify(aero_props(iaermod)%obj)
         514             :        end if
         515             :     end do
         516             :
         517           0 :     if (allocated(aero_props)) then
         518           0 :        deallocate(aero_props)
         519             :     endif
         520             :
         521           0 :   end subroutine aerosol_optics_cam_final
         522             :
         523             :   !===============================================================================
         524           0 :   subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, wa, ga, fa)
         525             :
         526             :     ! calculates aerosol sw radiative properties
         527             :
         528             :     integer,             intent(in) :: list_idx       ! index of the climate or a diagnostic list
         529             :     type(physics_state), intent(in), target :: state          ! state variables
         530             :
         531             :     type(physics_buffer_desc), pointer :: pbuf(:)
         532             :     integer,             intent(in) :: nnite          ! number of night columns
         533             :     integer,             intent(in) :: idxnite(nnite) ! local column indices of night columns
         534             :
         535             :     real(r8), intent(inout) :: tauxar(pcols,0:pver,nswbands) ! layer extinction optical depth
         536             :     real(r8), intent(inout) :: wa(pcols,0:pver,nswbands)     ! layer single-scatter albedo
         537             :     real(r8), intent(inout) :: ga(pcols,0:pver,nswbands)     ! asymmetry factor
         538             :     real(r8), intent(inout) :: fa(pcols,0:pver,nswbands)     ! forward scattered fraction
         539             :
         540             :     character(len=*), parameter :: prefix = 'aerosol_optics_cam_sw: '
         541             :
         542             :     integer :: ibin, nbins
         543             :     integer :: iwav, ilev
         544             :     integer :: icol, istat
         545             :     integer :: lchnk, ncol
         546             :
         547           0 :     type(aero_state_t), allocatable :: aero_state(:) ! array of aerosol state objects to allow for
         548             :                                                      ! multiple aerosol representations in the same sim
         549             :                                                      ! such as MAM and CARMA
         550             :
         551             :     class(aerosol_optics), pointer :: aero_optics
         552             :
         553             :     real(r8) :: dopaer(pcols)
         554             :     real(r8) :: mass(pcols,pver)
         555             :     real(r8) :: air_density(pcols,pver)
         556             :
         557           0 :     real(r8), allocatable :: pext(:)  ! parameterized specific extinction (m2/kg)
         558           0 :     real(r8), allocatable :: pabs(:)  ! parameterized specific absorption (m2/kg)
         559           0 :     real(r8), allocatable :: palb(:)  ! parameterized single scattering albedo
         560           0 :     real(r8), allocatable :: pasm(:)  ! parameterized asymmetry factor
         561             :
         562             :     character(len=ot_length) :: opticstype
         563             :     integer :: iaermod
         564             :
         565             :     real(r8) :: aodvis(pcols)              ! extinction optical depth in vis
         566             :     real(r8) :: aoduv(pcols)               ! extinction optical depth in uv
         567             :     real(r8) :: aodnir(pcols)              ! extinction optical depth in nir
         568             :     real(r8) :: absorb(pcols,pver)
         569             :     real(r8) :: aodabs(pcols)              ! absorption optical depth
         570             :
         571             :     real(r8) :: aodabsbc(pcols)             ! absorption optical depth of BC
         572             :
         573             :     real(r8) :: aodtot(pcols)
         574             :
         575             :     real(r8) :: extinct(pcols,pver)
         576             :     real(r8) :: extinctnir(pcols,pver)
         577             :     real(r8) :: extinctuv(pcols,pver)
         578             :
         579             :     real(r8) :: asymvis(pcols)              ! asymmetry factor * optical depth
         580             :     real(r8) :: asymext(pcols,pver)         ! asymmetry factor * extinction
         581             :
         582             :     real(r8) :: wetvol(pcols,pver)
         583             :     real(r8) :: watervol(pcols,pver)
         584             :
         585             :     real(r8) :: vol(pcols)
         586             :     real(r8) :: dustvol(pcols)
         587             :
         588             :     real(r8) :: scatdust(pcols)
         589             :     real(r8) :: absdust(pcols)
         590             :     real(r8) :: dustaodbin(pcols)
         591             :
         592             :     real(r8) :: scatbc(pcols)
         593             :     real(r8) :: absbc(pcols)
         594             :
         595             :     real(r8) :: scatpom(pcols)
         596             :     real(r8) :: abspom(pcols)
         597             :
         598             :     real(r8) :: scatsslt(pcols)
         599             :     real(r8) :: abssslt(pcols)
         600             :
         601             :     real(r8) :: scatsoa(pcols)
         602             :     real(r8) :: abssoa(pcols)
         603             :
         604             :     real(r8) :: scatsulf(pcols)
         605             :     real(r8) :: abssulf(pcols)
         606             :
         607             :     real(r8) :: burden(pcols)
         608             :     real(r8) :: burdendust(pcols), burdenso4(pcols), burdenbc(pcols), &
         609             :                 burdenpom(pcols), burdensoa(pcols), burdenseasalt(pcols)
         610             :
         611             :     real(r8) :: hygrodust(pcols), hygrosulf(pcols), hygrobc(pcols), &
         612             :                 hygropom(pcols), hygrosoa(pcols), hygrosslt(pcols)
         613             :
         614             :     real(r8) :: aodbin(pcols)
         615             :
         616           0 :     complex(r8), pointer :: specrefindex(:)     ! species refractive index
         617             :
         618             :     class(aerosol_state),      pointer :: aerostate
         619             :     class(aerosol_properties), pointer :: aeroprops
         620             :
         621             :     real(r8) :: specdens
         622             :     character(len=32) :: spectype            ! species type
         623           0 :     real(r8), pointer :: specmmr(:,:)
         624             :     real(r8)          :: hygro_aer           !
         625             :
         626             :     real(r8) :: scath2o, absh2o, sumscat, sumabs, sumhygro
         627             :
         628             :     real(r8) :: aodc                        ! aod of component
         629             :
         630             :     ! total species AOD
         631             :     real(r8) :: dustaod(pcols), sulfaod(pcols), bcaod(pcols), &
         632             :                 pomaod(pcols), soaaod(pcols), ssltaod(pcols)
         633             :
         634             :     real(r8) :: aodvisst(pcols) ! stratospheric extinction optical depth
         635             :     real(r8) :: aoduvst(pcols)  ! stratospheric extinction optical depth in uv
         636             :     real(r8) :: aodnirst(pcols) ! stratospheric extinction optical depth in nir
         637             :     real(r8) :: ssavis(pcols)
         638             :     integer :: troplev(pcols)
         639             :
         640           0 :     nullify(aero_optics)
         641             :
         642           0 :     lchnk = state%lchnk
         643           0 :     ncol  = state%ncol
         644             :
         645             :     !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         646           0 :     troplev(:) = 0
         647             :     !REMOVECAM_END
         648           0 :     call tropopause_findChemTrop(state, troplev)
         649             :
         650           0 :     mass(:ncol,:)        = state%pdeldry(:ncol,:)*rga
         651           0 :     air_density(:ncol,:) = state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
         652             :
         653           0 :     aodvis = 0._r8
         654           0 :     aodnir = 0._r8
         655           0 :     aoduv = 0._r8
         656           0 :     aodabs = 0._r8
         657           0 :     absorb = 0._r8
         658           0 :     aodtot = 0._r8
         659           0 :     tauxar = 0._r8
         660           0 :     extinct = 0._r8
         661           0 :     extinctnir = 0._r8
         662           0 :     extinctuv = 0._r8
         663           0 :     asymvis = 0.0_r8
         664           0 :     asymext = 0.0_r8
         665           0 :     ssavis = 0.0_r8
         666           0 :     aodvisst = 0.0_r8
         667           0 :     aoduvst = 0.0_r8
         668           0 :     aodnirst = 0.0_r8
         669             :
         670           0 :     burdendust = 0.0_r8
         671           0 :     burdenso4 = 0.0_r8
         672           0 :     burdenbc = 0.0_r8
         673           0 :     burdenpom = 0.0_r8
         674           0 :     burdensoa = 0.0_r8
         675           0 :     burdenseasalt = 0.0_r8
         676             :
         677           0 :     aodabsbc = 0.0_r8
         678           0 :     dustaod = 0.0_r8
         679           0 :     sulfaod = 0.0_r8
         680           0 :     pomaod = 0.0_r8
         681           0 :     soaaod = 0.0_r8
         682           0 :     bcaod = 0.0_r8
         683           0 :     ssltaod = 0.0_r8
         684             :
         685           0 :     if (num_aero_models<1) return
         686             :
         687           0 :     allocate(aero_state(num_aero_models), stat=istat)
         688           0 :     if (istat/=0) then
         689           0 :        call endrun(prefix//'array allocation error: aero_state')
         690             :     end if
         691             :
         692           0 :     iaermod = 0
         693           0 :     if (modal_active) then
         694           0 :        iaermod = iaermod+1
         695           0 :        aero_state(iaermod)%obj => modal_aerosol_state( state, pbuf )
         696             :     end if
         697             :
         698           0 :     allocate(pext(ncol), stat=istat)
         699           0 :     if (istat/=0) then
         700           0 :        call endrun(prefix//'array allocation error: pext')
         701             :     end if
         702           0 :     allocate(pabs(ncol), stat=istat)
         703           0 :     if (istat/=0) then
         704           0 :        call endrun(prefix//'array allocation error: pabs')
         705             :     end if
         706           0 :     allocate(palb(ncol), stat=istat)
         707           0 :     if (istat/=0) then
         708           0 :        call endrun(prefix//'array allocation error: palb')
         709             :     end if
         710           0 :     allocate(pasm(ncol), stat=istat)
         711           0 :     if (istat/=0) then
         712           0 :        call endrun(prefix//'array allocation error: pasm')
         713             :     end if
         714             :
         715           0 :     aeromodel: do iaermod = 1,num_aero_models
         716             :
         717           0 :        aeroprops => aero_props(iaermod)%obj
         718           0 :        aerostate => aero_state(iaermod)%obj
         719             :
         720           0 :        nbins=aeroprops%nbins(list_idx)
         721             :
         722           0 :        binloop: do ibin = 1, nbins
         723             :
         724           0 :           dustaodbin(:) = 0._r8
         725           0 :           burden(:) = 0._r8
         726           0 :           aodbin(:) = 0.0_r8
         727             :
         728           0 :           call aeroprops%optics_params(list_idx, ibin, opticstype=opticstype)
         729             :
         730           0 :           select case (trim(opticstype))
         731             :           case('modal') ! refractive method
         732             :              aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, &
         733           0 :                                                     ncol, pver, nswbands, nlwbands, crefwsw, crefwlw)
         734             :           case default
         735           0 :              call endrun(prefix//'optics method not recognized')
         736             :           end select
         737             :
         738           0 :           if (associated(aero_optics)) then
         739             :
         740           0 :              wetvol(:ncol,:pver) = aerostate%wet_volume(aeroprops, list_idx, ibin, ncol, pver)
         741           0 :              watervol(:ncol,:pver) = aerostate%water_volume(aeroprops, list_idx, ibin, ncol, pver)
         742             :
         743           0 :              wavelength: do iwav = 1, nswbands
         744             :
         745           0 :                 vertical: do ilev = 1, pver
         746             :
         747           0 :                    call aero_optics%sw_props(ncol, ilev, iwav, pext, pabs, palb, pasm )
         748             :
         749           0 :                    call init_diags
         750             :
         751           0 :                    column: do icol = 1,ncol
         752           0 :                       dopaer(icol) = pext(icol)*mass(icol,ilev)
         753           0 :                       tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol)
         754           0 :                       wa(icol,ilev,iwav) = wa(icol,ilev,iwav) + dopaer(icol)*palb(icol)
         755           0 :                       ga(icol,ilev,iwav) = ga(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol)
         756           0 :                       fa(icol,ilev,iwav) = fa(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol)*pasm(icol)
         757             :
         758           0 :                       call update_diags
         759             :
         760             :                    end do column
         761             :
         762             :                 end do vertical
         763             :              end do wavelength
         764             :
         765             :           else
         766           0 :              call endrun(prefix//'aero_optics object pointer not associated')
         767             :           end if
         768             :
         769           0 :           deallocate(aero_optics)
         770             :           nullify(aero_optics)
         771             :
         772           0 :           call output_bin_diags
         773             :
         774             :        end do binloop
         775             :     end do aeromodel
         776             :
         777           0 :     call output_tot_diags
         778             :
         779           0 :     deallocate(pext)
         780           0 :     deallocate(pabs)
         781           0 :     deallocate(palb)
         782           0 :     deallocate(pasm)
         783             :
         784           0 :     do iaermod = 1,num_aero_models
         785           0 :        deallocate(aero_state(iaermod)%obj)
         786           0 :        nullify(aero_state(iaermod)%obj)
         787             :     end do
         788             :
         789           0 :     deallocate(aero_state)
         790             :
         791             :   contains
         792             :
         793             :     !===============================================================================
         794           0 :     subroutine init_diags
         795           0 :       dustvol(:ncol)   = 0._r8
         796           0 :       scatdust(:ncol)  = 0._r8
         797           0 :       absdust(:ncol)   = 0._r8
         798           0 :       hygrodust(:ncol) = 0._r8
         799           0 :       scatsulf(:ncol)  = 0._r8
         800           0 :       abssulf(:ncol)   = 0._r8
         801           0 :       hygrosulf(:ncol) = 0._r8
         802           0 :       scatbc(:ncol)    = 0._r8
         803           0 :       absbc(:ncol)     = 0._r8
         804           0 :       hygrobc(:ncol)   = 0._r8
         805           0 :       scatpom(:ncol)   = 0._r8
         806           0 :       abspom(:ncol)    = 0._r8
         807           0 :       hygropom(:ncol)  = 0._r8
         808           0 :       scatsoa(:ncol)   = 0._r8
         809           0 :       abssoa(:ncol)    = 0._r8
         810           0 :       hygrosoa(:ncol)  = 0._r8
         811           0 :       scatsslt(:ncol)  = 0._r8
         812           0 :       abssslt(:ncol)   = 0._r8
         813           0 :       hygrosslt(:ncol) = 0._r8
         814           0 :     end subroutine init_diags
         815             :
         816             :     !===============================================================================
         817           0 :     subroutine update_diags
         818             :
         819             :       integer :: ispec
         820             :
         821           0 :       if (iwav==idx_uv_diag) then
         822           0 :          aoduv(icol) = aoduv(icol) + dopaer(icol)
         823           0 :          extinctuv(icol,ilev) = extinctuv(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
         824           0 :          if (ilev<=troplev(icol)) then
         825           0 :             aoduvst(icol) = aoduvst(icol) + dopaer(icol)
         826             :          end if
         827             :
         828           0 :       else if (iwav==idx_sw_diag) then ! vis
         829           0 :          aodvis(icol) = aodvis(icol) + dopaer(icol)
         830           0 :          aodabs(icol) = aodabs(icol) + pabs(icol)*mass(icol,ilev)
         831           0 :          extinct(icol,ilev) = extinct(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
         832           0 :          absorb(icol,ilev)  = absorb(icol,ilev) + pabs(icol)*air_density(icol,ilev)
         833           0 :          ssavis(icol)       = ssavis(icol) + dopaer(icol)*palb(icol)
         834           0 :          asymvis(icol)      = asymvis(icol) + dopaer(icol)*pasm(icol)
         835           0 :          asymext(icol,ilev) = asymext(icol,ilev) + dopaer(icol)*pasm(icol)*air_density(icol,ilev)/mass(icol,ilev)
         836             :
         837           0 :          aodbin(icol) = aodbin(icol) + dopaer(icol)
         838             :
         839           0 :          if (ilev<=troplev(icol)) then
         840           0 :             aodvisst(icol) = aodvisst(icol) + dopaer(icol)
         841             :          end if
         842             :
         843             :          ! loop over species ...
         844             :
         845           0 :          do ispec = 1, aeroprops%nspecies(list_idx,ibin)
         846             :             call aeroprops%get(ibin, ispec, list_ndx=list_idx, density=specdens, &
         847           0 :                  spectype=spectype, refindex_sw=specrefindex, hygro=hygro_aer)
         848           0 :             call aerostate%get_ambient_mmr(list_idx, ispec, ibin, specmmr)
         849             :
         850           0 :             burden(icol) = burden(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         851             :
         852           0 :             vol(icol) = specmmr(icol,ilev)/specdens
         853             :
         854           0 :             select case ( trim(spectype) )
         855             :             case('dust')
         856           0 :                dustvol(icol) = vol(icol)
         857           0 :                burdendust(icol) = burdendust(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         858           0 :                scatdust(icol) = vol(icol) * specrefindex(iwav)%re
         859           0 :                absdust(icol)  =-vol(icol) * specrefindex(iwav)%im
         860           0 :                hygrodust(icol)= vol(icol)*hygro_aer
         861             :             case('black-c')
         862           0 :                burdenbc(icol) = burdenbc(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         863           0 :                scatbc(icol) = vol(icol) * specrefindex(iwav)%re
         864           0 :                absbc(icol)  =-vol(icol) * specrefindex(iwav)%im
         865           0 :                hygrobc(icol)= vol(icol)*hygro_aer
         866             :             case('sulfate')
         867           0 :                burdenso4(icol) = burdenso4(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         868           0 :                scatsulf(icol) = vol(icol) * specrefindex(iwav)%re
         869           0 :                abssulf(icol)  =-vol(icol) * specrefindex(iwav)%im
         870           0 :                hygrosulf(icol)= vol(icol)*hygro_aer
         871             :             case('p-organic')
         872           0 :                burdenpom(icol) = burdenpom(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         873           0 :                scatpom(icol) = vol(icol) * specrefindex(iwav)%re
         874           0 :                abspom(icol)  =-vol(icol) * specrefindex(iwav)%im
         875           0 :                hygropom(icol)= vol(icol)*hygro_aer
         876             :             case('s-organic')
         877           0 :                burdensoa(icol) = burdensoa(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         878           0 :                scatsoa(icol) = vol(icol) * specrefindex(iwav)%re
         879           0 :                abssoa(icol) = -vol(icol) * specrefindex(iwav)%im
         880           0 :                hygrosoa(icol)= vol(icol)*hygro_aer
         881             :             case('seasalt')
         882           0 :                burdenseasalt(icol) = burdenseasalt(icol) + specmmr(icol,ilev)*mass(icol,ilev)
         883           0 :                scatsslt(icol) = vol(icol) * specrefindex(iwav)%re
         884           0 :                abssslt(icol) = -vol(icol) * specrefindex(iwav)%im
         885           0 :                hygrosslt(icol)= vol(icol)*hygro_aer
         886             :             end select
         887             :          end do
         888             :
         889           0 :          if (wetvol(icol,ilev)>1.e-40_r8 .and. vol(icol)>0._r8) then
         890             :
         891           0 :             dustaodbin(icol) = dustaodbin(icol) + dopaer(icol)*dustvol(icol)/wetvol(icol,ilev)
         892             :
         893             :             ! partition optical depth into contributions from each constituent
         894             :             ! assume contribution is proportional to refractive index X volume
         895             :
         896           0 :             scath2o = watervol(icol,ilev)*crefwsw(iwav)%re
         897           0 :             absh2o = -watervol(icol,ilev)*crefwsw(iwav)%im
         898             :             sumscat = scatsulf(icol) + scatpom(icol) + scatsoa(icol) + scatbc(icol) + &
         899           0 :                  scatdust(icol) + scatsslt(icol) + scath2o
         900             :             sumabs  = abssulf(icol) + abspom(icol) + abssoa(icol) + absbc(icol) + &
         901           0 :                  absdust(icol) + abssslt(icol) + absh2o
         902             :             sumhygro = hygrosulf(icol) + hygropom(icol) + hygrosoa(icol) + hygrobc(icol) + &
         903           0 :                  hygrodust(icol) + hygrosslt(icol)
         904             :
         905           0 :             scatdust(icol) = (scatdust(icol) + scath2o*hygrodust(icol)/sumhygro)/sumscat
         906           0 :             absdust(icol)  = (absdust(icol) + absh2o*hygrodust(icol)/sumhygro)/sumabs
         907             :
         908           0 :             scatsulf(icol) = (scatsulf(icol) + scath2o*hygrosulf(icol)/sumhygro)/sumscat
         909           0 :             abssulf(icol)  = (abssulf(icol) + absh2o*hygrosulf(icol)/sumhygro)/sumabs
         910             :
         911           0 :             scatpom(icol) = (scatpom(icol) + scath2o*hygropom(icol)/sumhygro)/sumscat
         912           0 :             abspom(icol)  = (abspom(icol) + absh2o*hygropom(icol)/sumhygro)/sumabs
         913             :
         914           0 :             scatsoa(icol) = (scatsoa(icol) + scath2o*hygrosoa(icol)/sumhygro)/sumscat
         915           0 :             abssoa(icol)  = (abssoa(icol) + absh2o*hygrosoa(icol)/sumhygro)/sumabs
         916             :
         917           0 :             scatbc(icol)= (scatbc(icol) + scath2o*hygrobc(icol)/sumhygro)/sumscat
         918           0 :             absbc(icol)  = (absbc(icol) +  absh2o*hygrobc(icol)/sumhygro)/sumabs
         919             :
         920           0 :             scatsslt(icol) = (scatsslt(icol) + scath2o*hygrosslt(icol)/sumhygro)/sumscat
         921           0 :             abssslt(icol)  = (abssslt(icol) + absh2o*hygrosslt(icol)/sumhygro)/sumabs
         922             :
         923             :
         924           0 :             aodabsbc(icol) = aodabsbc(icol) + absbc(icol)*dopaer(icol)*(1.0_r8-palb(icol))
         925             :
         926             :
         927             :
         928           0 :             aodc          = (absdust(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatdust(icol))*dopaer(icol)
         929           0 :             dustaod(icol) = dustaod(icol) + aodc
         930             :
         931           0 :             aodc          = (abssulf(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsulf(icol))*dopaer(icol)
         932           0 :             sulfaod(icol) = sulfaod(icol) + aodc
         933             :
         934           0 :             aodc          = (abspom(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatpom(icol))*dopaer(icol)
         935           0 :             pomaod(icol)  = pomaod(icol) + aodc
         936             :
         937           0 :             aodc          = (abssoa(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsoa(icol))*dopaer(icol)
         938           0 :             soaaod(icol)  = soaaod(icol) + aodc
         939             :
         940           0 :             aodc          = (absbc(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatbc(icol))*dopaer(icol)
         941           0 :             bcaod(icol)   = bcaod(icol) + aodc
         942             :
         943           0 :             aodc          = (abssslt(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsslt(icol))*dopaer(icol)
         944           0 :             ssltaod(icol) = ssltaod(icol) + aodc
         945             :
         946             :          end if
         947             :       else if (iwav==idx_nir_diag) then
         948             :          aodnir(icol) = aodnir(icol) + dopaer(icol)
         949             :          extinctnir(icol,ilev) = extinctnir(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
         950             :
         951             :          if (ilev<=troplev(icol)) then
         952             :             aodnirst(icol) = aodnirst(icol) + dopaer(icol)
         953             :          end if
         954             :
         955             :       end if
         956             :
         957           0 :       aodtot(icol) = aodtot(icol) + dopaer(icol)
         958             :
         959           0 :     end subroutine update_diags
         960             :
         961             :     !===============================================================================
         962           0 :     subroutine output_bin_diags
         963             :
         964             :       integer :: icol
         965             :
         966           0 :       if (list_idx == 0) then
         967             :
         968           0 :          call outfld(burdendn_fields(iaermod)%name(ibin), burden, pcols, lchnk)
         969           0 :          call outfld(aoddustdn_fields(iaermod)%name(ibin), dustaodbin, pcols, lchnk)
         970           0 :          call outfld(aodbindn_fields(iaermod)%name(ibin),      aodbin, pcols, lchnk)
         971             :
         972           0 :          do icol = 1, nnite
         973           0 :             burden(idxnite(icol))  = fillvalue
         974           0 :             aodbin(idxnite(icol)) = fillvalue
         975           0 :             dustaodbin(idxnite(icol)) = fillvalue
         976             :          end do
         977             :
         978           0 :          call outfld(burden_fields(iaermod)%name(ibin), burden, pcols, lchnk)
         979           0 :          call outfld(aoddust_fields(iaermod)%name(ibin), dustaodbin, pcols, lchnk)
         980           0 :          call outfld(aodbin_fields(iaermod)%name(ibin),      aodbin, pcols, lchnk)
         981             :
         982             :       endif
         983             :
         984           0 :     end subroutine output_bin_diags
         985             :
         986             :     !===============================================================================
         987           0 :     subroutine output_tot_diags
         988             :
         989             :       integer :: icol
         990             :
         991           0 :       call outfld('AODUVdn'//diag(list_idx),  aoduv,  pcols, lchnk)
         992           0 :       call outfld('AODVISdn'//diag(list_idx), aodvis, pcols, lchnk)
         993           0 :       call outfld('AODABSdn'//diag(list_idx),     aodabs,  pcols, lchnk)
         994           0 :       call outfld('AODNIRdn'//diag(list_idx), aodnir, pcols, lchnk)
         995           0 :       call outfld('AODTOTdn'//diag(list_idx), aodtot, pcols, lchnk)
         996           0 :       call outfld('EXTINCTUVdn'//diag(list_idx),  extinctuv,  pcols, lchnk)
         997           0 :       call outfld('EXTINCTNIRdn'//diag(list_idx), extinctnir, pcols, lchnk)
         998           0 :       call outfld('EXTINCTdn'//diag(list_idx),  extinct,  pcols, lchnk)
         999           0 :       call outfld('ABSORBdn'//diag(list_idx),   absorb,  pcols, lchnk)
        1000           0 :       call outfld('EXTxASYMdn'//diag(list_idx), asymext, pcols, lchnk)
        1001           0 :       call outfld('AODxASYMdn'//diag(list_idx), asymvis, pcols, lchnk)
        1002             :
        1003           0 :       call outfld('AODVISstdn'//diag(list_idx), aodvisst,pcols, lchnk)
        1004           0 :       call outfld('AODUVstdn'//diag(list_idx),  aoduvst, pcols, lchnk)
        1005           0 :       call outfld('AODNIRstdn'//diag(list_idx), aodnirst,pcols, lchnk)
        1006             :
        1007           0 :       do icol = 1, nnite
        1008           0 :          aodvis(idxnite(icol)) = fillvalue
        1009           0 :          aodnir(idxnite(icol)) = fillvalue
        1010           0 :          aoduv(idxnite(icol)) = fillvalue
        1011           0 :          aodabs(idxnite(icol)) = fillvalue
        1012           0 :          aodtot(idxnite(icol)) = fillvalue
        1013           0 :          extinct(idxnite(icol),:) = fillvalue
        1014           0 :          extinctnir(idxnite(icol),:) = fillvalue
        1015           0 :          extinctuv(idxnite(icol),:) = fillvalue
        1016           0 :          absorb(idxnite(icol),:)  = fillvalue
        1017           0 :          asymext(idxnite(icol),:) = fillvalue
        1018           0 :          asymvis(idxnite(icol)) = fillvalue
        1019           0 :          aodabs(idxnite(icol))    = fillvalue
        1020           0 :          aodvisst(idxnite(icol))  = fillvalue
        1021           0 :          aoduvst(idxnite(icol))   = fillvalue
        1022           0 :          aodnirst(idxnite(icol))  = fillvalue
        1023             :       end do
        1024             :
        1025           0 :       call outfld('AODUV'//diag(list_idx),  aoduv,  pcols, lchnk)
        1026           0 :       call outfld('AODVIS'//diag(list_idx), aodvis, pcols, lchnk)
        1027           0 :       call outfld('AODABS'//diag(list_idx), aodabs,  pcols, lchnk)
        1028           0 :       call outfld('AODNIR'//diag(list_idx), aodnir, pcols, lchnk)
        1029           0 :       call outfld('AODTOT'//diag(list_idx), aodtot, pcols, lchnk)
        1030           0 :       call outfld('EXTINCTUV'//diag(list_idx),  extinctuv,  pcols, lchnk)
        1031           0 :       call outfld('EXTINCTNIR'//diag(list_idx), extinctnir, pcols, lchnk)
        1032           0 :       call outfld('EXTINCT'//diag(list_idx),  extinct,  pcols, lchnk)
        1033           0 :       call outfld('ABSORB'//diag(list_idx),   absorb,  pcols, lchnk)
        1034           0 :       call outfld('EXTxASYM'//diag(list_idx), asymext, pcols, lchnk)
        1035           0 :       call outfld('AODxASYM'//diag(list_idx), asymvis, pcols, lchnk)
        1036           0 :       call outfld('AODVISst'//diag(list_idx), aodvisst,pcols, lchnk)
        1037           0 :       call outfld('AODUVst'//diag(list_idx),  aoduvst, pcols, lchnk)
        1038           0 :       call outfld('AODNIRst'//diag(list_idx), aodnirst,pcols, lchnk)
        1039             :
        1040             :       ! These diagnostics are output only for climate list
        1041           0 :       if (list_idx == 0) then
        1042           0 :          do icol = 1, ncol
        1043           0 :             if (aodvis(icol) > 1.e-10_r8) then
        1044           0 :                ssavis(icol) = ssavis(icol)/aodvis(icol)
        1045             :             else
        1046           0 :                ssavis(icol) = 0.925_r8
        1047             :             endif
        1048             :          end do
        1049           0 :          call outfld('SSAVISdn',        ssavis,        pcols, lchnk)
        1050             :
        1051           0 :          call outfld('BURDENDUSTdn',    burdendust,    pcols, lchnk)
        1052           0 :          call outfld('BURDENSO4dn' ,    burdenso4,     pcols, lchnk)
        1053           0 :          call outfld('BURDENPOMdn' ,    burdenpom,     pcols, lchnk)
        1054           0 :          call outfld('BURDENSOAdn' ,    burdensoa,     pcols, lchnk)
        1055           0 :          call outfld('BURDENBCdn'  ,    burdenbc,      pcols, lchnk)
        1056           0 :          call outfld('BURDENSEASALTdn', burdenseasalt, pcols, lchnk)
        1057             :
        1058           0 :          call outfld('AODABSBCdn',      aodabsbc,      pcols, lchnk)
        1059             :
        1060           0 :          call outfld('AODDUSTdn',       dustaod,       pcols, lchnk)
        1061           0 :          call outfld('AODSO4dn',        sulfaod,       pcols, lchnk)
        1062           0 :          call outfld('AODPOMdn',        pomaod,        pcols, lchnk)
        1063           0 :          call outfld('AODSOAdn',        soaaod,        pcols, lchnk)
        1064           0 :          call outfld('AODBCdn',         bcaod,         pcols, lchnk)
        1065           0 :          call outfld('AODSSdn',         ssltaod,       pcols, lchnk)
        1066             :
        1067             :
        1068           0 :          do icol = 1, nnite
        1069             :
        1070           0 :             ssavis(idxnite(icol))     = fillvalue
        1071           0 :             asymvis(idxnite(icol))    = fillvalue
        1072             :
        1073           0 :             burdendust(idxnite(icol)) = fillvalue
        1074           0 :             burdenso4(idxnite(icol))  = fillvalue
        1075           0 :             burdenpom(idxnite(icol))  = fillvalue
        1076           0 :             burdensoa(idxnite(icol))  = fillvalue
        1077           0 :             burdenbc(idxnite(icol))   = fillvalue
        1078           0 :             burdenseasalt(idxnite(icol)) = fillvalue
        1079           0 :             aodabsbc(idxnite(icol)) = fillvalue
        1080             :
        1081           0 :             dustaod(idxnite(icol))    = fillvalue
        1082           0 :             sulfaod(idxnite(icol))     = fillvalue
        1083           0 :             pomaod(idxnite(icol))     = fillvalue
        1084           0 :             soaaod(idxnite(icol))     = fillvalue
        1085           0 :             bcaod(idxnite(icol))      = fillvalue
        1086           0 :             ssltaod(idxnite(icol)) = fillvalue
        1087             :
        1088             :          end do
        1089             :
        1090           0 :          call outfld('SSAVIS',        ssavis,        pcols, lchnk)
        1091           0 :          call outfld('AODxASYM',      asymvis,       pcols, lchnk)
        1092           0 :          call outfld('BURDENDUST',    burdendust,    pcols, lchnk)
        1093           0 :          call outfld('BURDENSO4' ,    burdenso4,     pcols, lchnk)
        1094           0 :          call outfld('BURDENPOM' ,    burdenpom,     pcols, lchnk)
        1095           0 :          call outfld('BURDENSOA' ,    burdensoa,     pcols, lchnk)
        1096           0 :          call outfld('BURDENBC'  ,    burdenbc,      pcols, lchnk)
        1097           0 :          call outfld('BURDENSEASALT', burdenseasalt, pcols, lchnk)
        1098           0 :          call outfld('AODABSBC',      aodabsbc,      pcols, lchnk)
        1099           0 :          call outfld('AODDUST',       dustaod,       pcols, lchnk)
        1100           0 :          call outfld('AODSO4',        sulfaod,       pcols, lchnk)
        1101           0 :          call outfld('AODPOM',        pomaod,        pcols, lchnk)
        1102           0 :          call outfld('AODSOA',        soaaod,        pcols, lchnk)
        1103           0 :          call outfld('AODBC',         bcaod,         pcols, lchnk)
        1104           0 :          call outfld('AODSS',         ssltaod,       pcols, lchnk)
        1105             :
        1106             :       end if
        1107             :
        1108           0 :     end subroutine output_tot_diags
        1109             :
        1110             :   end subroutine aerosol_optics_cam_sw
        1111             :
        1112             :   !===============================================================================
        1113           0 :   subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar)
        1114             :
        1115             :     ! calculates aerosol lw radiative properties
        1116             :
        1117             :     integer,             intent(in)  :: list_idx ! index of the climate or a diagnostic list
        1118             :     type(physics_state), intent(in), target :: state    ! state variables
        1119             :
        1120             :     type(physics_buffer_desc), pointer :: pbuf(:)
        1121             :
        1122             :     real(r8), intent(inout) :: tauxar(pcols,pver,nlwbands) ! layer absorption optical depth
        1123             :
        1124             :
        1125             :     real(r8) :: dopaer(pcols)
        1126             :     real(r8) :: mass(pcols,pver)
        1127             :
        1128             :     character(len=*), parameter :: prefix = 'aerosol_optics_cam_lw: '
        1129             :
        1130             :     integer :: ibin, nbins
        1131             :     integer :: iwav, ilev
        1132             :     integer :: ncol, icol, istat
        1133             :
        1134           0 :     type(aero_state_t), allocatable :: aero_state(:) ! array of aerosol state objects to allow for
        1135             :                                                      ! multiple aerosol representations in the same sim
        1136             :                                                      ! such as MAM and CARMA
        1137             :
        1138             :     class(aerosol_optics), pointer :: aero_optics
        1139             :     class(aerosol_state),      pointer :: aerostate
        1140             :     class(aerosol_properties), pointer :: aeroprops
        1141             :
        1142           0 :     real(r8), allocatable :: pabs(:)
        1143             :
        1144             :     character(len=32) :: opticstype
        1145             :     integer :: iaermod
        1146             :
        1147             :     real(r8) :: lwabs(pcols,pver)
        1148           0 :     lwabs = 0._r8
        1149           0 :     tauxar = 0._r8
        1150             :
        1151           0 :     nullify(aero_optics)
        1152             :
        1153           0 :     allocate(aero_state(num_aero_models), stat=istat)
        1154           0 :     if (istat/=0) then
        1155           0 :        call endrun(prefix//'array allocation error: aero_state')
        1156             :     end if
        1157             :
        1158           0 :     iaermod = 0
        1159           0 :     if (modal_active) then
        1160           0 :        iaermod = iaermod+1
        1161           0 :        aero_state(iaermod)%obj => modal_aerosol_state( state, pbuf )
        1162             :     end if
        1163             :
        1164           0 :     ncol = state%ncol
        1165             :
        1166           0 :     mass(:ncol,:) = state%pdeldry(:ncol,:)*rga
        1167             :
        1168           0 :     allocate(pabs(ncol), stat=istat)
        1169           0 :     if (istat/=0) then
        1170           0 :        call endrun(prefix//'array allocation error: pabs')
        1171             :     end if
        1172             :
        1173           0 :     aeromodel: do iaermod = 1,num_aero_models
        1174             :
        1175           0 :        aeroprops => aero_props(iaermod)%obj
        1176           0 :        aerostate => aero_state(iaermod)%obj
        1177             :
        1178           0 :        nbins=aero_props(iaermod)%obj%nbins(list_idx)
        1179             :
        1180           0 :        binloop: do ibin = 1, nbins
        1181             :
        1182           0 :           call aeroprops%optics_params(list_idx, ibin, opticstype=opticstype)
        1183             :
        1184           0 :           select case (trim(opticstype))
        1185             :           case('modal') ! refractive method
        1186             :              aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, &
        1187           0 :                                                     ncol, pver, nswbands, nlwbands, crefwsw, crefwlw)
        1188             :           case default
        1189           0 :              call endrun(prefix//'optics method not recognized')
        1190             :           end select
        1191             :
        1192           0 :           if (associated(aero_optics)) then
        1193             :
        1194           0 :              wavelength: do iwav = 1, nlwbands
        1195             :
        1196           0 :                 vertical: do ilev = 1, pver
        1197           0 :                    call aero_optics%lw_props(ncol, ilev, iwav, pabs )
        1198             :
        1199           0 :                    column: do icol = 1, ncol
        1200           0 :                       dopaer(icol) = pabs(icol)*mass(icol,ilev)
        1201           0 :                       tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol)
        1202           0 :                       lwabs(icol,ilev) = lwabs(icol,ilev) + pabs(icol)
        1203             :                    end do column
        1204             :
        1205             :                 end do vertical
        1206             :
        1207             :              end do wavelength
        1208             :
        1209             :           else
        1210           0 :              call endrun(prefix//'aero_optics object pointer not associated')
        1211             :           end if
        1212             :
        1213           0 :           deallocate(aero_optics)
        1214           0 :           nullify(aero_optics)
        1215             :
        1216             :        end do binloop
        1217             :     end do aeromodel
        1218             :
        1219           0 :     call outfld('TOTABSLW'//diag(list_idx),  lwabs(:,:), pcols, state%lchnk)
        1220             :
        1221           0 :     if (lw10um_indx>0) then
        1222           0 :        call outfld('AODABSLW'//diag(list_idx), tauxar(:,:,lw10um_indx), pcols, state%lchnk)
        1223             :     end if
        1224             :
        1225           0 :     deallocate(pabs)
        1226             :
        1227           0 :     do iaermod = 1,num_aero_models
        1228           0 :        deallocate(aero_state(iaermod)%obj)
        1229           0 :        nullify(aero_state(iaermod)%obj)
        1230             :     end do
        1231             :
        1232           0 :     deallocate(aero_state)
        1233             :
        1234           0 :   end subroutine aerosol_optics_cam_lw
        1235             :
        1236             :   !===============================================================================
        1237             :   ! Private routines
        1238             :   !===============================================================================
        1239             :
        1240           0 :   subroutine read_water_refindex(infilename)
        1241             :     use cam_pio_utils, only: cam_pio_openfile
        1242             :     use pio, only: file_desc_t, var_desc_t, pio_inq_dimlen, pio_inq_dimid, pio_inq_varid, &
        1243             :                    pio_get_var, PIO_NOWRITE, pio_closefile, pio_noerr
        1244             :
        1245             :
        1246             :     ! read water refractive index file and set module data
        1247             :
        1248             :     character*(*), intent(in) :: infilename   ! modal optics filename
        1249             :
        1250             :     ! Local variables
        1251             :
        1252             :     integer            :: i, ierr
        1253             :     type(file_desc_t)  :: ncid              ! pio file handle
        1254             :     integer            :: did               ! dimension ids
        1255             :     integer            :: dimlen            ! dimension lengths
        1256             :     type(var_desc_t)   :: vid               ! variable ids
        1257             :     real(r8) :: refrwsw(nswbands), refiwsw(nswbands) ! real, imaginary ref index for water visible
        1258             :     real(r8) :: refrwlw(nlwbands), refiwlw(nlwbands) ! real, imaginary ref index for water infrared
        1259             :
        1260             :     character(len=*), parameter :: prefix = 'read_water_refindex: '
        1261             :     !----------------------------------------------------------------------------
        1262             :
        1263             :     ! open file
        1264           0 :     call cam_pio_openfile(ncid, infilename, PIO_NOWRITE)
        1265             :
        1266             :     ! inquire dimensions.  Check that file values match parameter values.
        1267             :
        1268           0 :     ierr = pio_inq_dimid(ncid, 'lw_band', did)
        1269           0 :     if (ierr /= pio_noerr ) then
        1270           0 :        call endrun(prefix//'pio_inq_dimid lw_band')
        1271             :     end if
        1272           0 :     ierr = pio_inq_dimlen(ncid, did, dimlen)
        1273           0 :     if (ierr /= pio_noerr ) then
        1274           0 :        call endrun(prefix//'pio_inq_dimlen lw_band')
        1275             :     end if
        1276           0 :     if (dimlen /= nlwbands) then
        1277           0 :        write(iulog,*) 'lw_band len=', dimlen, ' from ', infilename, ' ne nlwbands=', nlwbands
        1278           0 :        call endrun(prefix//'bad lw_band value')
        1279             :     endif
        1280             :
        1281           0 :     ierr = pio_inq_dimid(ncid, 'sw_band', did)
        1282           0 :     if (ierr /= pio_noerr ) then
        1283           0 :        call endrun(prefix//'pio_inq_dimid sw_band')
        1284             :     end if
        1285           0 :     ierr = pio_inq_dimlen(ncid, did, dimlen)
        1286           0 :     if (ierr /= pio_noerr ) then
        1287           0 :        call endrun(prefix//'pio_inq_dimlen sw_band')
        1288             :     end if
        1289           0 :     if (dimlen /= nswbands) then
        1290           0 :        write(iulog,*) 'sw_band len=', dimlen, ' from ', infilename, ' ne nswbands=', nswbands
        1291           0 :        call endrun(prefix//'bad sw_band value')
        1292             :     endif
        1293             :
        1294             :     ! read variables
        1295           0 :     ierr = pio_inq_varid(ncid, 'refindex_real_water_sw', vid)
        1296           0 :     if (ierr /= pio_noerr ) then
        1297           0 :        call endrun(prefix//'pio_inq_varid refindex_real_water_sw')
        1298             :     end if
        1299           0 :     ierr = pio_get_var(ncid, vid, refrwsw)
        1300           0 :     if (ierr /= pio_noerr ) then
        1301           0 :        call endrun(prefix//'pio_get_var refrwsw')
        1302             :     end if
        1303             :
        1304           0 :     ierr = pio_inq_varid(ncid, 'refindex_im_water_sw', vid)
        1305           0 :     if (ierr /= pio_noerr ) then
        1306           0 :        call endrun(prefix//'pio_inq_varid refindex_im_water_sw')
        1307             :     end if
        1308           0 :     ierr = pio_get_var(ncid, vid, refiwsw)
        1309           0 :     if (ierr /= pio_noerr ) then
        1310           0 :        call endrun(prefix//'pio_get_var refiwsw')
        1311             :     end if
        1312             :
        1313           0 :     ierr = pio_inq_varid(ncid, 'refindex_real_water_lw', vid)
        1314           0 :     if (ierr /= pio_noerr ) then
        1315           0 :        call endrun(prefix//'pio_inq_varid refindex_real_water_lw')
        1316             :     end if
        1317           0 :     ierr = pio_get_var(ncid, vid, refrwlw)
        1318           0 :     if (ierr /= pio_noerr ) then
        1319           0 :        call endrun(prefix//'pio_get_var refrwlw')
        1320             :     end if
        1321             :
        1322           0 :     ierr = pio_inq_varid(ncid, 'refindex_im_water_lw', vid)
        1323           0 :     if (ierr /= pio_noerr ) then
        1324           0 :        call endrun(prefix//'pio_inq_varid refindex_im_water_lw')
        1325             :     end if
        1326           0 :     ierr = pio_get_var(ncid, vid, refiwlw)
        1327           0 :     if (ierr /= pio_noerr ) then
        1328           0 :        call endrun(prefix//'pio_get_var refiwlw')
        1329             :     end if
        1330             :
        1331             :     ! set complex representation of refractive indices as module data
        1332           0 :     do i = 1, nswbands
        1333           0 :        crefwsw(i) = cmplx(refrwsw(i), abs(refiwsw(i)), kind=r8)
        1334             :     end do
        1335           0 :     do i = 1, nlwbands
        1336           0 :        crefwlw(i) = cmplx(refrwlw(i), abs(refiwlw(i)), kind=r8)
        1337             :     end do
        1338             :
        1339           0 :     call pio_closefile(ncid)
        1340             :
        1341           0 :   end subroutine read_water_refindex
        1342             :
        1343           0 : end module aerosol_optics_cam

Generated by: LCOV version 1.14