LCOV - code coverage report
Current view: | top level - physics/cam - microp_aero.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 49 | 308 | 15.9 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 3 | 8 | 37.5 %

              Line data    Source code
           1             : module microp_aero
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !   CAM driver layer for aerosol activation processes.
           6             : !
           7             : ! ***N.B.*** This module is currently hardcoded to recognize only the aerosols/modes that
           8             : !            affect the climate calculation.  This is implemented by using list
           9             : !            index 0 in all the calls to rad_constituent interfaces.
          10             : !
          11             : ! Author: Andrew Gettelman
          12             : ! Based on code from: Hugh Morrison, Xiaohong Liu and Steve Ghan
          13             : ! May 2010
          14             : ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
          15             : !                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
          16             : ! for questions contact Andrew Gettelman  (andrew@ucar.edu)
          17             : ! Modifications: A. Gettelman Nov 2010  - changed to support separation of
          18             : !                  microphysics and macrophysics and concentrate aerosol information here
          19             : !                B. Eaton, Sep 2014 - Refactored to move CAM interface code into the CAM
          20             : !                  interface modules and preserve just the driver layer functionality here.
          21             : !
          22             : !---------------------------------------------------------------------------------
          23             :
          24             : use shr_kind_mod,     only: r8=>shr_kind_r8
          25             : use spmd_utils,       only: masterproc
          26             : use ppgrid,           only: pcols, pver, pverp, begchunk, endchunk
          27             : use ref_pres,         only: top_lev => trop_cloud_top_lev
          28             : use physconst,        only: rair
          29             : use constituents,     only: cnst_get_ind
          30             : use physics_types,    only: physics_state, physics_ptend, physics_ptend_init, physics_ptend_sum, &
          31             :                             physics_state_copy, physics_update
          32             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field, &
          33             :                             pbuf_get_chunk
          34             : use phys_control,     only: phys_getopts, use_hetfrz_classnuc
          35             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
          36             :                             rad_cnst_get_mode_num
          37             :
          38             : use nucleate_ice_cam, only: use_preexisting_ice, nucleate_ice_cam_readnl, nucleate_ice_cam_register, &
          39             :                             nucleate_ice_cam_init, nucleate_ice_cam_calc
          40             :
          41             : use ndrop,            only: ndrop_init, dropmixnuc
          42             : use ndrop_bam,        only: ndrop_bam_init, ndrop_bam_run, ndrop_bam_ccn
          43             :
          44             : use hetfrz_classnuc_cam, only: hetfrz_classnuc_cam_readnl, hetfrz_classnuc_cam_register, hetfrz_classnuc_cam_init, &
          45             :                                hetfrz_classnuc_cam_calc
          46             :
          47             : use cam_history,      only: addfld, add_default, outfld
          48             : use cam_logfile,      only: iulog
          49             : use cam_abortutils,       only: endrun
          50             :
          51             : use aerosol_properties_mod, only: aerosol_properties
          52             : use modal_aerosol_properties_mod, only: modal_aerosol_properties
          53             :
          54             : use aerosol_state_mod, only: aerosol_state
          55             : use modal_aerosol_state_mod, only: modal_aerosol_state
          56             :
          57             : implicit none
          58             : private
          59             : save
          60             :
          61             : public :: microp_aero_init, microp_aero_run, microp_aero_readnl, microp_aero_register
          62             : public :: microp_aero_final
          63             : public :: aerosol_state_object
          64             : public :: aerosol_properties_object
          65             :
          66             : ! Private module data
          67             : character(len=16)   :: eddy_scheme
          68             : real(r8), parameter :: unset_r8 = huge(1.0_r8)
          69             :
          70             : ! contact freezing due to dust
          71             : ! dust number mean radius (m), Zender et al JGR 2003 assuming number mode radius of 0.6 micron, sigma=2
          72             : real(r8), parameter :: rn_dst1 = 0.258e-6_r8
          73             : real(r8), parameter :: rn_dst2 = 0.717e-6_r8
          74             : real(r8), parameter :: rn_dst3 = 1.576e-6_r8
          75             : real(r8), parameter :: rn_dst4 = 3.026e-6_r8
          76             :
          77             : ! Namelist parameters
          78             : real(r8) :: bulk_scale    ! prescribed aerosol bulk sulfur scale factor
          79             : real(r8) :: npccn_scale   ! scaling for activated number
          80             : real(r8) :: wsub_scale    ! scaling for sub-grid vertical velocity (liquid)
          81             : real(r8) :: wsubi_scale   ! scaling for sub-grid vertical velocity (ice)
          82             : real(r8) :: wsub_min      ! minimum sub-grid vertical velocity (liquid) before scale factor
          83             : real(r8) :: wsub_min_asf  ! minimum sub-grid vertical velocity (liquid) after scale factor
          84             : real(r8) :: wsubi_min     ! minimum sub-grid vertical velocity (ice)
          85             :
          86             : ! smallest mixing ratio considered in microphysics
          87             : real(r8), parameter :: qsmall = 1.e-18_r8
          88             :
          89             : ! minimum allowed cloud fraction
          90             : real(r8), parameter :: mincld = 0.0001_r8
          91             :
          92             : ! indices in state%q and pbuf structures
          93             : integer :: cldliq_idx = -1
          94             : integer :: cldice_idx = -1
          95             : integer :: numliq_idx = -1
          96             : integer :: numice_idx = -1
          97             : integer :: kvh_idx = -1
          98             : integer :: tke_idx = -1
          99             : integer :: wp2_idx = -1
         100             : integer :: ast_idx = -1
         101             : integer :: cldo_idx = -1
         102             : integer :: dgnumwet_idx = -1
         103             :
         104             : ! Bulk aerosols
         105             : character(len=20), allocatable :: aername(:)
         106             : real(r8), allocatable :: num_to_mass_aer(:)
         107             :
         108             : integer :: naer_all      ! number of aerosols affecting climate
         109             : integer :: idxsul   = -1 ! index in aerosol list for sulfate
         110             : integer :: idxdst2  = -1 ! index in aerosol list for dust2
         111             : integer :: idxdst3  = -1 ! index in aerosol list for dust3
         112             : integer :: idxdst4  = -1 ! index in aerosol list for dust4
         113             :
         114             : ! modal aerosols
         115             : logical :: clim_modal_aero
         116             :
         117             : integer :: mode_accum_idx  = -1  ! index of accumulation mode
         118             : integer :: mode_aitken_idx = -1  ! index of aitken mode
         119             : integer :: mode_coarse_idx = -1  ! index of coarse mode
         120             : integer :: mode_coarse_dst_idx = -1  ! index of coarse dust mode
         121             : integer :: mode_coarse_slt_idx = -1  ! index of coarse sea salt mode
         122             : integer :: coarse_dust_idx = -1  ! index of dust in coarse mode
         123             : integer :: coarse_nacl_idx = -1  ! index of nacl in coarse mode
         124             : integer :: coarse_so4_idx = -1  ! index of sulfate in coarse mode
         125             :
         126             : integer :: npccn_idx, rndst_idx, nacon_idx
         127             :
         128             : logical  :: separate_dust = .false.
         129             :
         130             : type aero_state_t
         131             :    class(aerosol_state), pointer :: obj=>null()
         132             : end type aero_state_t
         133             :
         134             : class(aerosol_properties), pointer :: aero_props_obj=>null()
         135             : type(aero_state_t), pointer :: aero_state(:) => null()
         136             :
         137             : !=========================================================================================
         138             : contains
         139             : !=========================================================================================
         140             :
         141           0 : subroutine microp_aero_register
         142             :    !-----------------------------------------------------------------------
         143             :    !
         144             :    ! Purpose:
         145             :    ! Register pbuf fields for aerosols needed by microphysics
         146             :    !
         147             :    ! Author: Cheryl Craig October 2012
         148             :    !
         149             :    !-----------------------------------------------------------------------
         150             :    use ppgrid,         only: pcols
         151             :    use physics_buffer, only: pbuf_add_field, dtype_r8
         152             :
         153           0 :    call pbuf_add_field('NPCCN',      'physpkg',dtype_r8,(/pcols,pver/), npccn_idx)
         154             :
         155           0 :    call pbuf_add_field('RNDST',      'physpkg',dtype_r8,(/pcols,pver,4/), rndst_idx)
         156           0 :    call pbuf_add_field('NACON',      'physpkg',dtype_r8,(/pcols,pver,4/), nacon_idx)
         157             :
         158           0 :    call nucleate_ice_cam_register()
         159           0 :    call hetfrz_classnuc_cam_register()
         160             :
         161           0 : end subroutine microp_aero_register
         162             :
         163             : !=========================================================================================
         164             :
         165           0 : subroutine microp_aero_init(phys_state,pbuf2d)
         166             :
         167             :    !-----------------------------------------------------------------------
         168             :    !
         169             :    ! Purpose:
         170             :    ! Initialize constants for aerosols needed by microphysics
         171             :    !
         172             :    ! Author: Andrew Gettelman May 2010
         173             :    !
         174             :    !-----------------------------------------------------------------------
         175             :
         176             :    type(physics_state), pointer       :: phys_state(:)
         177             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         178             :
         179             :    ! local variables
         180             :    integer  :: iaer, ierr
         181             :    integer  :: m, n, nmodes, nspec
         182             :
         183             :    character(len=32) :: str32
         184             :    character(len=*), parameter :: routine = 'microp_aero_init'
         185             :    logical :: history_amwg
         186           0 :    type(physics_buffer_desc), pointer :: pbuf(:)
         187             :    integer :: c
         188             :
         189             :    !-----------------------------------------------------------------------
         190             :
         191             :    ! Query the PBL eddy scheme
         192             :    call phys_getopts(eddy_scheme_out          = eddy_scheme,  &
         193           0 :                      history_amwg_out         = history_amwg )
         194             :
         195             :    ! Access the physical properties of the aerosols that are affecting the climate
         196             :    ! by using routines from the rad_constituents module.
         197             :
         198             :    ! get indices into state and pbuf structures
         199           0 :    call cnst_get_ind('CLDLIQ', cldliq_idx)
         200           0 :    call cnst_get_ind('CLDICE', cldice_idx)
         201           0 :    call cnst_get_ind('NUMLIQ', numliq_idx)
         202           0 :    call cnst_get_ind('NUMICE', numice_idx)
         203             :
         204           0 :    select case(trim(eddy_scheme))
         205             :    case ('diag_TKE')
         206           0 :       tke_idx      = pbuf_get_index('tke')
         207             :    case ('CLUBB_SGS')
         208           0 :       wp2_idx = pbuf_get_index('WP2_nadv')
         209             :    case default
         210           0 :       kvh_idx      = pbuf_get_index('kvh')
         211             :    end select
         212             :
         213             :    ! clim_modal_aero determines whether modal aerosols are used in the climate calculation.
         214             :    ! The modal aerosols can be either prognostic or prescribed.
         215           0 :    call rad_cnst_get_info(0, nmodes=nmodes)
         216           0 :    clim_modal_aero = (nmodes > 0)
         217             :
         218           0 :    ast_idx = pbuf_get_index('AST')
         219             :
         220           0 :    if (clim_modal_aero) then
         221             :
         222           0 :       cldo_idx = pbuf_get_index('CLDO')
         223           0 :       dgnumwet_idx = pbuf_get_index('DGNUMWET')
         224             :
         225           0 :       aero_props_obj => modal_aerosol_properties()
         226           0 :       if (.not.associated(aero_props_obj)) then
         227           0 :          call endrun('ma_convproc_init: construction of modal_aerosol_properties object failed')
         228             :       end if
         229           0 :       call ndrop_init(aero_props_obj)
         230           0 :       call nucleate_ice_cam_init(mincld, bulk_scale, pbuf2d, aero_props=aero_props_obj)
         231             :
         232           0 :       allocate(aero_state(begchunk:endchunk))
         233           0 :       do c = begchunk,endchunk
         234           0 :          pbuf => pbuf_get_chunk(pbuf2d, c)
         235           0 :          aero_state(c)%obj => modal_aerosol_state( phys_state(c), pbuf )
         236           0 :          if (.not.associated(aero_state(c)%obj)) then
         237           0 :             call endrun('microp_aero_init: construction of modal_aerosol_state object failed')
         238             :          end if
         239             :       end do
         240             :
         241             :       ! Init indices for specific modes/species
         242             :
         243             :       ! mode index for specified mode types
         244           0 :       do m = 1, nmodes
         245           0 :          call rad_cnst_get_info(0, m, mode_type=str32)
         246           0 :          select case (trim(str32))
         247             :          case ('accum')
         248           0 :             mode_accum_idx = m
         249             :          case ('aitken')
         250           0 :             mode_aitken_idx = m
         251             :          case ('coarse')
         252           0 :             mode_coarse_idx = m
         253             :          case ('coarse_dust')
         254           0 :             mode_coarse_dst_idx = m
         255             :          case ('coarse_seasalt')
         256           0 :             mode_coarse_slt_idx = m
         257             :          end select
         258             :       end do
         259             :
         260             :       ! check if coarse dust is in separate mode
         261           0 :       separate_dust = mode_coarse_dst_idx > 0
         262             :
         263             :       ! for 3-mode
         264           0 :       if ( mode_coarse_dst_idx<0 ) mode_coarse_dst_idx = mode_coarse_idx
         265           0 :       if ( mode_coarse_slt_idx<0 ) mode_coarse_slt_idx = mode_coarse_idx
         266             :
         267             :       ! Check that required mode types were found
         268             :       if (mode_accum_idx == -1 .or. mode_aitken_idx == -1 .or. &
         269           0 :           mode_coarse_dst_idx == -1.or. mode_coarse_slt_idx == -1) then
         270           0 :          write(iulog,*) routine//': ERROR required mode type not found - mode idx:', &
         271           0 :             mode_accum_idx, mode_aitken_idx, mode_coarse_dst_idx, mode_coarse_slt_idx
         272           0 :          call endrun(routine//': ERROR required mode type not found')
         273             :       end if
         274             :
         275             :       ! species indices for specified types
         276             :       ! find indices for the dust and seasalt species in the coarse mode
         277           0 :       call rad_cnst_get_info(0, mode_coarse_dst_idx, nspec=nspec)
         278           0 :       do n = 1, nspec
         279           0 :          call rad_cnst_get_info(0, mode_coarse_dst_idx, n, spec_type=str32)
         280           0 :          select case (trim(str32))
         281             :          case ('dust')
         282           0 :             coarse_dust_idx = n
         283             :          end select
         284             :       end do
         285           0 :       call rad_cnst_get_info(0, mode_coarse_slt_idx, nspec=nspec)
         286           0 :       do n = 1, nspec
         287           0 :          call rad_cnst_get_info(0, mode_coarse_slt_idx, n, spec_type=str32)
         288           0 :          select case (trim(str32))
         289             :          case ('seasalt')
         290           0 :             coarse_nacl_idx = n
         291             :          end select
         292             :       end do
         293           0 :       if (mode_coarse_idx>0) then
         294           0 :          call rad_cnst_get_info(0, mode_coarse_idx, nspec=nspec)
         295           0 :          do n = 1, nspec
         296           0 :             call rad_cnst_get_info(0, mode_coarse_idx, n, spec_type=str32)
         297           0 :             select case (trim(str32))
         298             :             case ('sulfate')
         299           0 :                coarse_so4_idx = n
         300             :             end select
         301             :          end do
         302             :       endif
         303             :
         304             :       ! Check that required mode specie types were found
         305           0 :       if ( coarse_dust_idx == -1 .or. coarse_nacl_idx == -1 ) then
         306           0 :          write(iulog,*) routine//': ERROR required mode-species type not found - indicies:', &
         307           0 :             coarse_dust_idx, coarse_nacl_idx
         308           0 :          call endrun(routine//': ERROR required mode-species type not found')
         309             :       end if
         310             :
         311             :    else
         312             :
         313             :       ! Props needed for BAM number concentration calcs.
         314             :
         315           0 :       call rad_cnst_get_info(0, naero=naer_all)
         316             :       allocate( &
         317           0 :          aername(naer_all),        &
         318           0 :          num_to_mass_aer(naer_all) )
         319             :
         320           0 :       do iaer = 1, naer_all
         321             :          call rad_cnst_get_aer_props(0, iaer, &
         322           0 :             aername         = aername(iaer), &
         323           0 :             num_to_mass_aer = num_to_mass_aer(iaer) )
         324             :
         325             :          ! Look for sulfate, dust, and soot in this list (Bulk aerosol only)
         326           0 :          if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
         327           0 :          if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
         328           0 :          if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
         329           0 :          if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
         330             :       end do
         331             :
         332           0 :       call ndrop_bam_init()
         333           0 :       call nucleate_ice_cam_init(mincld, bulk_scale, pbuf2d)
         334             :
         335             :    end if
         336             :
         337           0 :    call addfld('LCLOUD', (/ 'lev' /), 'A', ' ',   'Liquid cloud fraction used in stratus activation')
         338             :
         339           0 :    call addfld('WSUB',   (/ 'lev' /), 'A', 'm/s', 'Diagnostic sub-grid vertical velocity'                   )
         340           0 :    call addfld('WSUBI',  (/ 'lev' /), 'A', 'm/s', 'Diagnostic sub-grid vertical velocity for ice'           )
         341             :
         342           0 :    if (history_amwg) then
         343           0 :       call add_default ('WSUB     ', 1, ' ')
         344             :    end if
         345             :
         346           0 :    if (use_hetfrz_classnuc) then
         347           0 :       if (associated(aero_props_obj)) then
         348           0 :          call hetfrz_classnuc_cam_init(mincld, aero_props_obj)
         349             :       else
         350           0 :          call endrun(routine//': cannot use hetfrz_classnuc without prognostic aerosols')
         351             :       endif
         352             :    endif
         353             :
         354           0 : end subroutine microp_aero_init
         355             :
         356             : !=========================================================================================
         357             : ! returns a pointer to an aerosol state object for a given chunk index
         358           0 : function aerosol_state_object(lchnk) result(obj)
         359             :
         360             :   integer,intent(in) :: lchnk ! local chunk index
         361             :   class(aerosol_state), pointer :: obj ! aerosol state object pointer for local chunk
         362             :
         363           0 :   obj => aero_state(lchnk)%obj
         364             :
         365           0 : end function aerosol_state_object
         366             :
         367             : !=========================================================================================
         368             : ! returns a pointer to an aerosol properties object
         369        1536 : function aerosol_properties_object() result(obj)
         370             :
         371             :   class(aerosol_properties), pointer :: obj ! aerosol properties object pointer
         372             :
         373        1536 :   obj => aero_props_obj
         374             :
         375        1536 : end function aerosol_properties_object
         376             :
         377             : !=========================================================================================
         378             :
         379        1536 : subroutine microp_aero_final
         380             :
         381             :   integer :: c
         382             :
         383        1536 :   if (associated(aero_props_obj)) then
         384           0 :      deallocate(aero_props_obj)
         385             :   end if
         386        1536 :   nullify(aero_props_obj)
         387             :
         388        1536 :   if (associated(aero_state)) then
         389           0 :      do c = begchunk,endchunk
         390           0 :         deallocate(aero_state(c)%obj)
         391             :      end do
         392           0 :      deallocate(aero_state)
         393             :      nullify(aero_state)
         394             :   end if
         395             :
         396        1536 : end subroutine microp_aero_final
         397             :
         398             : !=========================================================================================
         399             :
         400        1536 : subroutine microp_aero_readnl(nlfile)
         401             :
         402             :    use namelist_utils,  only: find_group_name
         403             :    use units,           only: getunit, freeunit
         404             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_real8
         405             :
         406             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         407             :
         408             :    ! Namelist variables
         409             :    real(r8) :: microp_aero_bulk_scale = unset_r8   ! prescribed aerosol bulk sulfur scale factor
         410             :    real(r8) :: microp_aero_npccn_scale = unset_r8  ! prescribed aerosol bulk sulfur scale factor
         411             :    real(r8) :: microp_aero_wsub_scale = unset_r8   ! subgrid vertical velocity (liquid) scale factor
         412             :    real(r8) :: microp_aero_wsubi_scale = unset_r8  ! subgrid vertical velocity (ice) scale factor
         413             :    real(r8) :: microp_aero_wsub_min = unset_r8     ! subgrid vertical velocity (liquid) minimum (before scale factor)
         414             :    real(r8) :: microp_aero_wsub_min_asf = unset_r8 ! subgrid vertical velocity (liquid) minimum (after scale factor)
         415             :    real(r8) :: microp_aero_wsubi_min = unset_r8    ! subgrid vertical velocity (ice) minimum
         416             :
         417             :    ! Local variables
         418             :    integer :: unitn, ierr
         419             :    character(len=*), parameter :: subname = 'microp_aero_readnl'
         420             :
         421             :    namelist /microp_aero_nl/ microp_aero_bulk_scale, microp_aero_npccn_scale, microp_aero_wsub_min, &
         422             :                              microp_aero_wsubi_min, microp_aero_wsub_scale, microp_aero_wsubi_scale, microp_aero_wsub_min_asf
         423             :    !-----------------------------------------------------------------------------
         424             :
         425        1536 :    if (masterproc) then
         426           2 :       unitn = getunit()
         427           2 :       open( unitn, file=trim(nlfile), status='old' )
         428           2 :       call find_group_name(unitn, 'microp_aero_nl', status=ierr)
         429           2 :       if (ierr == 0) then
         430           2 :          read(unitn, microp_aero_nl, iostat=ierr)
         431           2 :          if (ierr /= 0) then
         432           0 :             call endrun(subname // ':: ERROR reading namelist')
         433             :          end if
         434             :       end if
         435           2 :       close(unitn)
         436           2 :       call freeunit(unitn)
         437             :    end if
         438             :
         439             :    ! Broadcast namelist variables
         440        1536 :    call mpi_bcast(microp_aero_bulk_scale, 1, mpi_real8, mstrid, mpicom, ierr)
         441        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_bulk_scale")
         442        1536 :    call mpi_bcast(microp_aero_npccn_scale, 1, mpi_real8, mstrid, mpicom, ierr)
         443        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_npccn_scale")
         444        1536 :    call mpi_bcast(microp_aero_wsub_scale, 1, mpi_real8, mstrid, mpicom, ierr)
         445        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_scale")
         446        1536 :    call mpi_bcast(microp_aero_wsubi_scale, 1, mpi_real8, mstrid, mpicom, ierr)
         447        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsubi_scale")
         448        1536 :    call mpi_bcast(microp_aero_wsub_min, 1, mpi_real8, mstrid, mpicom, ierr)
         449        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_min")
         450        1536 :    call mpi_bcast(microp_aero_wsub_min_asf, 1, mpi_real8, mstrid, mpicom, ierr)
         451        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsub_min_asf")
         452        1536 :    call mpi_bcast(microp_aero_wsubi_min, 1, mpi_real8, mstrid, mpicom, ierr)
         453        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: microp_aero_wsubi_min")
         454             :
         455             :    ! set local variables
         456        1536 :    bulk_scale = microp_aero_bulk_scale
         457        1536 :    npccn_scale = microp_aero_npccn_scale
         458        1536 :    wsub_scale = microp_aero_wsub_scale
         459        1536 :    wsubi_scale = microp_aero_wsubi_scale
         460        1536 :    wsub_min = microp_aero_wsub_min
         461        1536 :    wsub_min_asf = microp_aero_wsub_min_asf
         462        1536 :    wsubi_min = microp_aero_wsubi_min
         463             :
         464        1536 :    if(bulk_scale == unset_r8) call endrun(subname//": FATAL: bulk_scale is not set")
         465        1536 :    if(npccn_scale == unset_r8) call endrun(subname//": FATAL: npccn_scale is not set")
         466        1536 :    if(wsub_scale == unset_r8) call endrun(subname//": FATAL: wsub_scale is not set")
         467        1536 :    if(wsubi_scale == unset_r8) call endrun(subname//": FATAL: wsubi_scale is not set")
         468        1536 :    if(wsub_min == unset_r8) call endrun(subname//": FATAL: wsub_min is not set")
         469        1536 :    if(wsub_min_asf == unset_r8) call endrun(subname//": FATAL: wsub_min_asf is not set")
         470        1536 :    if(wsubi_min == unset_r8) call endrun(subname//": FATAL: wsubi_min is not set")
         471             :
         472        1536 :    call nucleate_ice_cam_readnl(nlfile)
         473        1536 :    call hetfrz_classnuc_cam_readnl(nlfile)
         474             :
         475        1536 : end subroutine microp_aero_readnl
         476             :
         477             : !=========================================================================================
         478             :
         479           0 : subroutine microp_aero_run ( &
         480             :    state, ptend_all, deltatin, pbuf)
         481             :
         482             :    ! input arguments
         483             :    type(physics_state),         intent(in)    :: state
         484             :    type(physics_ptend),         intent(out)   :: ptend_all
         485             :    real(r8),                    intent(in)    :: deltatin     ! time step (s)
         486             :    type(physics_buffer_desc),   pointer       :: pbuf(:)
         487             :
         488             :    ! local workspace
         489             :    ! all units mks unless otherwise stated
         490             :
         491             :    integer :: i, k, m
         492             :    integer :: itim_old
         493             :
         494           0 :    type(physics_state), target :: state1                ! Local copy of state variable
         495           0 :    type(physics_ptend) :: ptend_loc
         496             :
         497           0 :    real(r8), pointer :: ast(:,:)
         498             :
         499           0 :    real(r8), pointer :: npccn(:,:)      ! number of CCN (liquid activated)
         500             :
         501           0 :    real(r8), pointer :: rndst(:,:,:)    ! radius of 4 dust bins for contact freezing
         502           0 :    real(r8), pointer :: nacon(:,:,:)    ! number in 4 dust bins for contact freezing
         503             :
         504           0 :    real(r8), pointer :: num_coarse(:,:) ! number m.r. of coarse mode
         505           0 :    real(r8), pointer :: coarse_dust(:,:) ! mass m.r. of coarse dust
         506           0 :    real(r8), pointer :: coarse_nacl(:,:) ! mass m.r. of coarse nacl
         507           0 :    real(r8), pointer :: coarse_so4(:,:)  ! mass m.r. of coarse sulfate
         508             :
         509           0 :    real(r8), pointer :: kvh(:,:)        ! vertical eddy diff coef (m2 s-1)
         510           0 :    real(r8), pointer :: tke(:,:)        ! TKE from the UW PBL scheme (m2 s-2)
         511           0 :    real(r8), pointer :: wp2(:,:)        ! CLUBB vertical velocity variance
         512             :
         513           0 :    real(r8), pointer :: cldn(:,:)       ! cloud fraction
         514           0 :    real(r8), pointer :: cldo(:,:)       ! old cloud fraction
         515             :
         516           0 :    real(r8), pointer :: dgnumwet(:,:,:) ! aerosol mode diameter
         517             :
         518           0 :    real(r8), pointer :: aer_mmr(:,:)    ! aerosol mass mixing ratio
         519             :
         520             :    real(r8) :: rho(pcols,pver)     ! air density (kg m-3)
         521             :
         522             :    real(r8) :: lcldm(pcols,pver)   ! liq cloud fraction
         523             :
         524             :    real(r8) :: lcldn(pcols,pver)   ! fractional coverage of new liquid cloud
         525             :    real(r8) :: lcldo(pcols,pver)   ! fractional coverage of old liquid cloud
         526             :    real(r8) :: cldliqf(pcols,pver) ! fractional of total cloud that is liquid
         527             :    real(r8) :: qcld                ! total cloud water
         528             :    real(r8) :: nctend_mixnuc(pcols,pver)
         529             :    real(r8) :: dum, dum2           ! temporary dummy variable
         530             :    real(r8) :: dmc, ssmc, so4mc    ! variables for modal scheme.
         531             :
         532             :    ! bulk aerosol variables
         533           0 :    real(r8), allocatable :: naer2(:,:,:)    ! bulk aerosol number concentration (1/m3)
         534           0 :    real(r8), allocatable :: maerosol(:,:,:) ! bulk aerosol mass conc (kg/m3)
         535             :
         536             :    real(r8) :: wsub(pcols,pver)    ! diagnosed sub-grid vertical velocity st. dev. (m/s)
         537             :    real(r8) :: wsubi(pcols,pver)   ! diagnosed sub-grid vertical velocity ice (m/s)
         538             :    real(r8) :: nucboas
         539             :
         540             :    real(r8) :: wght
         541             :
         542             :    integer :: lchnk, ncol, astat
         543             :
         544           0 :    real(r8), allocatable :: factnum(:,:,:) ! activation fraction for aerosol number
         545             :
         546             :    class(aerosol_state), pointer :: aero_state1_obj
         547             :
         548             :    !-------------------------------------------------------------------------------
         549             :
         550           0 :    nullify(aero_state1_obj)
         551             :
         552           0 :    call physics_state_copy(state,state1)
         553             :
         554           0 :    lchnk = state1%lchnk
         555           0 :    ncol  = state1%ncol
         556             :
         557           0 :    itim_old = pbuf_old_tim_idx()
         558           0 :    call pbuf_get_field(pbuf, ast_idx,      ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         559             :
         560           0 :    call pbuf_get_field(pbuf, npccn_idx, npccn)
         561             :
         562           0 :    call pbuf_get_field(pbuf, nacon_idx, nacon)
         563           0 :    call pbuf_get_field(pbuf, rndst_idx, rndst)
         564             :
         565           0 :    call physics_ptend_init(ptend_all, state%psetcols, 'microp_aero')
         566             :
         567           0 :    if (clim_modal_aero) then
         568             :       ! create an aerosol state object specifically for cam state1
         569           0 :       aero_state1_obj => modal_aerosol_state( state1, pbuf )
         570           0 :       if (.not.associated(aero_state1_obj)) then
         571           0 :          call endrun('microp_aero_run: construction of aero_state1_obj modal_aerosol_state object failed')
         572             :       end if
         573             :
         574           0 :       itim_old = pbuf_old_tim_idx()
         575             :
         576           0 :       call pbuf_get_field(pbuf, ast_idx,  cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         577           0 :       call pbuf_get_field(pbuf, cldo_idx, cldo, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
         578           0 :       call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet)
         579             :    end if
         580             :
         581             :    ! initialize output
         582           0 :    npccn(1:ncol,1:pver)    = 0._r8
         583             :
         584           0 :    nacon(1:ncol,1:pver,:)  = 0._r8
         585             :
         586             :    ! set default or fixed dust bins for contact freezing
         587           0 :    rndst(1:ncol,1:pver,1) = rn_dst1
         588           0 :    rndst(1:ncol,1:pver,2) = rn_dst2
         589           0 :    rndst(1:ncol,1:pver,3) = rn_dst3
         590           0 :    rndst(1:ncol,1:pver,4) = rn_dst4
         591             :
         592             :    ! initialize time-varying parameters
         593           0 :    do k = top_lev, pver
         594           0 :       do i = 1, ncol
         595           0 :          rho(i,k) = state1%pmid(i,k)/(rair*state1%t(i,k))
         596             :       end do
         597             :    end do
         598             :
         599           0 :    if (clim_modal_aero) then
         600             :       ! mode number mixing ratios
         601           0 :       call rad_cnst_get_mode_num(0, mode_coarse_dst_idx, 'a', state1, pbuf, num_coarse)
         602             :
         603             :       ! mode specie mass m.r.
         604           0 :       call rad_cnst_get_aer_mmr(0, mode_coarse_dst_idx, coarse_dust_idx, 'a', state1, pbuf, coarse_dust)
         605           0 :       call rad_cnst_get_aer_mmr(0, mode_coarse_slt_idx, coarse_nacl_idx, 'a', state1, pbuf, coarse_nacl)
         606           0 :       if (mode_coarse_idx>0) then
         607           0 :          call rad_cnst_get_aer_mmr(0, mode_coarse_idx, coarse_so4_idx, 'a', state1, pbuf, coarse_so4)
         608             :       endif
         609             :
         610             :    else
         611             :       ! init number/mass arrays for bulk aerosols
         612             :       allocate( &
         613           0 :          naer2(pcols,pver,naer_all), &
         614           0 :          maerosol(pcols,pver,naer_all))
         615             :
         616           0 :       do m = 1, naer_all
         617           0 :          call rad_cnst_get_aer_mmr(0, m, state1, pbuf, aer_mmr)
         618           0 :          maerosol(:ncol,:,m) = aer_mmr(:ncol,:)*rho(:ncol,:)
         619             :
         620           0 :          if (m .eq. idxsul) then
         621           0 :             naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)*bulk_scale
         622             :          else
         623           0 :             naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)
         624             :          end if
         625             :       end do
         626             :    end if
         627             :
         628             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         629             :    ! More refined computation of sub-grid vertical velocity
         630             :    ! Set to be zero at the surface by initialization.
         631             :
         632           0 :    select case (trim(eddy_scheme))
         633             :    case ('diag_TKE')
         634           0 :       call pbuf_get_field(pbuf, tke_idx, tke)
         635             :    case ('CLUBB_SGS')
         636           0 :       itim_old = pbuf_old_tim_idx()
         637           0 :       call pbuf_get_field(pbuf, wp2_idx, wp2, start=(/1,1,itim_old/),kount=(/pcols,pverp,1/))
         638           0 :       allocate(tke(pcols,pverp))
         639           0 :       tke(:ncol,:) = (3._r8/2._r8)*wp2(:ncol,:)
         640             :
         641             :    case default
         642           0 :       call pbuf_get_field(pbuf, kvh_idx, kvh)
         643             :    end select
         644             :
         645             :    ! Set minimum values above top_lev.
         646           0 :    wsub(:ncol,:top_lev-1)  = wsub_min
         647           0 :    wsubi(:ncol,:top_lev-1) = wsubi_min
         648             :
         649           0 :    do k = top_lev, pver
         650           0 :       do i = 1, ncol
         651             :
         652           0 :          select case (trim(eddy_scheme))
         653             :          case ('diag_TKE', 'CLUBB_SGS')
         654           0 :             wsub(i,k) = sqrt(0.5_r8*(tke(i,k) + tke(i,k+1))*(2._r8/3._r8))
         655           0 :             wsub(i,k) = min(wsub(i,k),10._r8)
         656             :          case default
         657             :             ! get sub-grid vertical velocity from diff coef.
         658             :             ! following morrison et al. 2005, JAS
         659             :             ! assume mixing length of 30 m
         660           0 :             dum = (kvh(i,k) + kvh(i,k+1))/2._r8/30._r8
         661             :             ! use maximum sub-grid vertical vel of 10 m/s
         662           0 :             dum = min(dum, 10._r8)
         663             :             ! set wsub to value at current vertical level
         664           0 :             wsub(i,k)  = dum
         665             :          end select
         666             :
         667           0 :          wsubi(i,k) = max(wsubi_min, wsub(i,k)) * wsubi_scale
         668           0 :          if (.not. use_preexisting_ice) then
         669           0 :             wsubi(i,k) = min(wsubi(i,k), 0.2_r8)
         670             :          endif
         671             :
         672           0 :          wsub(i,k)  = max(wsub_min, wsub(i,k)) * wsub_scale
         673             :
         674             :       end do
         675             :    end do
         676             :
         677           0 :    call outfld('WSUB',   wsub, pcols, lchnk)
         678           0 :    call outfld('WSUBI', wsubi, pcols, lchnk)
         679             :
         680           0 :    if (trim(eddy_scheme) == 'CLUBB_SGS') deallocate(tke)
         681             :
         682             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         683             :    !ICE Nucleation
         684             :
         685           0 :    if (associated(aero_props_obj).and.associated(aero_state1_obj)) then
         686           0 :       call nucleate_ice_cam_calc(state1, wsubi, pbuf, deltatin, ptend_loc, aero_props_obj, aero_state1_obj)
         687             :    else
         688           0 :       call nucleate_ice_cam_calc(state1, wsubi, pbuf, deltatin, ptend_loc)
         689             :    end if
         690             :
         691           0 :    call physics_ptend_sum(ptend_loc, ptend_all, ncol)
         692           0 :    call physics_update(state1, ptend_loc, deltatin)
         693             :
         694             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         695             :    ! get liquid cloud fraction, check for minimum
         696             :
         697           0 :    do k = top_lev, pver
         698           0 :       do i = 1, ncol
         699           0 :          lcldm(i,k) = max(ast(i,k), mincld)
         700             :       end do
         701             :    end do
         702             :
         703             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         704             :    ! Droplet Activation
         705             :
         706           0 :    if (clim_modal_aero) then
         707             :
         708             :       ! for modal aerosol
         709             :
         710             :       ! partition cloud fraction into liquid water part
         711           0 :       lcldn = 0._r8
         712           0 :       lcldo = 0._r8
         713           0 :       cldliqf = 0._r8
         714           0 :       do k = top_lev, pver
         715           0 :          do i = 1, ncol
         716           0 :             qcld = state1%q(i,k,cldliq_idx) + state1%q(i,k,cldice_idx)
         717           0 :             if (qcld > qsmall) then
         718           0 :                lcldn(i,k)   = cldn(i,k)*state1%q(i,k,cldliq_idx)/qcld
         719           0 :                lcldo(i,k)   = cldo(i,k)*state1%q(i,k,cldliq_idx)/qcld
         720           0 :                cldliqf(i,k) = state1%q(i,k,cldliq_idx)/qcld
         721             :             end if
         722             :          end do
         723             :       end do
         724             :
         725           0 :       call outfld('LCLOUD', lcldn, pcols, lchnk)
         726             :
         727           0 :       allocate(factnum(pcols,pver,aero_props_obj%nbins()),stat=astat)
         728           0 :       if (astat/=0) then
         729           0 :          call endrun('microp_aero_run: not able to allocate factnum')
         730             :       endif
         731             :
         732             :       ! If not using preexsiting ice, then only use cloudbourne aerosol for the
         733             :       ! liquid clouds. This is the same behavior as CAM5.
         734           0 :       if (use_preexisting_ice) then
         735             :          call dropmixnuc( aero_props_obj, aero_state1_obj, &
         736             :               state1, ptend_loc, deltatin, pbuf, wsub, wsub_min_asf, &
         737           0 :               cldn, cldo, cldliqf, nctend_mixnuc, factnum)
         738             :       else
         739           0 :          cldliqf = 1._r8
         740             :          call dropmixnuc( aero_props_obj, aero_state1_obj, &
         741             :               state1, ptend_loc, deltatin, pbuf, wsub, wsub_min_asf, &
         742           0 :               lcldn, lcldo, cldliqf, nctend_mixnuc, factnum)
         743             :       end if
         744             :
         745           0 :       npccn(:ncol,:) = nctend_mixnuc(:ncol,:)
         746             :
         747           0 :       npccn(:ncol,:) = npccn(:ncol,:) * npccn_scale
         748             :
         749             :    else
         750             :
         751             :       ! for bulk aerosol
         752             :
         753             :       ! no tendencies returned from ndrop_bam_run, so just init ptend here
         754           0 :       call physics_ptend_init(ptend_loc, state1%psetcols, 'none')
         755             :
         756           0 :       do k = top_lev, pver
         757           0 :          do i = 1, ncol
         758             :
         759           0 :             if (state1%q(i,k,cldliq_idx) >= qsmall) then
         760             :
         761             :                ! get droplet activation rate
         762             :
         763             :                call ndrop_bam_run( &
         764           0 :                   wsub(i,k), state1%t(i,k), rho(i,k), naer2(i,k,:), naer_all, &
         765           0 :                   naer_all, maerosol(i,k,:),  &
         766           0 :                   dum2)
         767           0 :                dum = dum2
         768             :             else
         769             :                dum = 0._r8
         770             :             end if
         771             :
         772           0 :             npccn(i,k) = (dum*lcldm(i,k) - state1%q(i,k,numliq_idx))/deltatin
         773             :          end do
         774             :       end do
         775             :
         776             :    end if
         777             :
         778           0 :    call physics_ptend_sum(ptend_loc, ptend_all, ncol)
         779           0 :    call physics_update(state1, ptend_loc, deltatin)
         780             :
         781             :
         782             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         783             :    ! Contact freezing  (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
         784             :    ! estimate rndst and nanco for 4 dust bins here to pass to MG microphysics
         785             :
         786           0 :    do k = top_lev, pver
         787           0 :       do i = 1, ncol
         788             :
         789           0 :          if (state1%t(i,k) < 269.15_r8) then
         790             :
         791           0 :             if (clim_modal_aero) then
         792             :
         793             :                ! For modal aerosols:
         794             :                !  use size '3' for dust coarse mode...
         795             :                !  scale by dust fraction in coarse mode
         796             :
         797           0 :                dmc  = coarse_dust(i,k)
         798           0 :                ssmc = coarse_nacl(i,k)
         799             :
         800           0 :                if ( separate_dust ) then
         801             :                   ! 7-mode -- has separate dust and seasalt mode types and no need for weighting
         802             :                   wght = 1._r8
         803             :                else
         804           0 :                   so4mc = coarse_so4(i,k)
         805             :                   ! 3-mode -- needs weighting for dust since dust, seasalt, and sulfate  are combined in the "coarse" mode type
         806           0 :                   wght = dmc/(ssmc + dmc + so4mc)
         807             :                endif
         808             :
         809           0 :                if (dmc > 0.0_r8) then
         810           0 :                   nacon(i,k,3) = wght*num_coarse(i,k)*rho(i,k)
         811             :                else
         812           0 :                   nacon(i,k,3) = 0._r8
         813             :                end if
         814             :
         815             :                !also redefine parameters based on size...
         816             :
         817           0 :                rndst(i,k,3) = 0.5_r8*dgnumwet(i,k,mode_coarse_dst_idx)
         818           0 :                if (rndst(i,k,3) <= 0._r8) then
         819           0 :                   rndst(i,k,3) = rn_dst3
         820             :                end if
         821             :
         822             :             else
         823             :
         824             :                !For Bulk Aerosols: set equal to aerosol number for dust for bins 2-4 (bin 1=0)
         825             :
         826           0 :                if (idxdst2 > 0) then
         827           0 :                   nacon(i,k,2) = naer2(i,k,idxdst2)
         828             :                end if
         829           0 :                if (idxdst3 > 0) then
         830           0 :                   nacon(i,k,3) = naer2(i,k,idxdst3)
         831             :                end if
         832           0 :                if (idxdst4 > 0) then
         833           0 :                   nacon(i,k,4) = naer2(i,k,idxdst4)
         834             :                end if
         835             :             end if
         836             :
         837             :          end if
         838             :       end do
         839             :    end do
         840             :
         841             :    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         842             :    !bulk aerosol ccn concentration (modal does it in ndrop, from dropmixnuc)
         843             :
         844           0 :    if (.not. clim_modal_aero) then
         845             :
         846             :       ! ccn concentration as diagnostic
         847           0 :       call ndrop_bam_ccn(lchnk, ncol, maerosol, naer2)
         848             :
         849           0 :       deallocate( &
         850             :          naer2,    &
         851           0 :          maerosol)
         852             :
         853             :    end if
         854             :
         855             :    ! heterogeneous freezing
         856           0 :    if (use_hetfrz_classnuc) then
         857             :
         858           0 :       call hetfrz_classnuc_cam_calc(aero_props_obj, aero_state1_obj, state1, deltatin, factnum, pbuf)
         859             :
         860             :    end if
         861             :
         862           0 :    if (clim_modal_aero) then
         863           0 :       deallocate(factnum)
         864             :    end if
         865             :
         866           0 :    if (associated(aero_state1_obj)) then
         867             :       ! destroy the aerosol state object
         868           0 :       deallocate(aero_state1_obj)
         869             :       nullify(aero_state1_obj)
         870             :    endif
         871             :
         872           0 :  end subroutine microp_aero_run
         873             :
         874             : !=========================================================================================
         875             :
         876           0 : end module microp_aero

Generated by: LCOV version 1.14