LCOV - code coverage report
Current view: | top level - physics/cam - nucleate_ice_cam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 18 | 383 | 4.7 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 4 | 25.0 %

              Line data    Source code
           1             : module nucleate_ice_cam
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : !
           5             : !  CAM Interfaces for nucleate_ice module.
           6             : !
           7             : !  B. Eaton - Sept 2014
           8             : !---------------------------------------------------------------------------------
           9             :
          10             : use shr_kind_mod,   only: r8=>shr_kind_r8
          11             : use spmd_utils,     only: masterproc
          12             : use ppgrid,         only: pcols, pver
          13             : use physconst,      only: pi, rair, tmelt
          14             : use constituents,   only: pcnst, cnst_get_ind
          15             : use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
          16             : use physics_buffer, only: physics_buffer_desc
          17             : use phys_control,   only: use_hetfrz_classnuc
          18             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props
          19             :
          20             : use physics_buffer, only: pbuf_add_field, dtype_r8, pbuf_old_tim_idx, &
          21             :                           pbuf_get_index, pbuf_get_field, &
          22             :                           pbuf_set_field
          23             : use cam_history,    only: addfld, add_default, outfld
          24             :
          25             : use ref_pres,       only: top_lev => trop_cloud_top_lev
          26             : use wv_saturation,  only: qsat_water
          27             :
          28             : use cam_logfile,    only: iulog
          29             : use cam_abortutils, only: endrun
          30             :
          31             : use nucleate_ice,   only: nucleati_init, nucleati
          32             :
          33             : use aerosol_properties_mod, only: aerosol_properties
          34             : use aerosol_state_mod, only: aerosol_state
          35             :
          36             : use phys_control,   only: cam_physpkg_is
          37             :
          38             : implicit none
          39             : private
          40             : save
          41             :
          42             : public :: &
          43             :    nucleate_ice_cam_readnl,   &
          44             :    nucleate_ice_cam_register, &
          45             :    nucleate_ice_cam_init,     &
          46             :    nucleate_ice_cam_calc
          47             :
          48             : ! Namelist variables
          49             : logical, public, protected :: use_preexisting_ice = .false.
          50             : logical                    :: hist_preexisting_ice = .false.
          51             : logical                    :: nucleate_ice_incloud = .false.
          52             : logical                    :: nucleate_ice_use_troplev = .false.
          53             : real(r8)                   :: nucleate_ice_subgrid = -1._r8
          54             : real(r8)                   :: nucleate_ice_subgrid_strat = -1._r8
          55             : real(r8)                   :: nucleate_ice_strat = 0.0_r8
          56             :
          57             : ! Vars set via init method.
          58             : real(r8) :: mincld      ! minimum allowed cloud fraction
          59             : real(r8) :: bulk_scale  ! prescribed aerosol bulk sulfur scale factor
          60             :
          61             : ! constituent indices
          62             : integer :: &
          63             :    cldliq_idx = -1, &
          64             :    cldice_idx = -1, &
          65             :    numice_idx = -1
          66             :
          67             : integer :: &
          68             :    naai_idx = -1,     &
          69             :    naai_hom_idx = -1
          70             :
          71             : integer :: &
          72             :    ast_idx   = -1
          73             :
          74             : integer :: &
          75             :     qsatfac_idx = -1
          76             :
          77             : ! Bulk aerosols
          78             : character(len=20), allocatable :: aername(:)
          79             : real(r8), allocatable :: num_to_mass_aer(:)
          80             :
          81             : integer :: naer_all = -1 ! number of aerosols affecting climate
          82             : integer :: idxsul   = -1 ! index in aerosol list for sulfate
          83             : integer :: idxdst1  = -1 ! index in aerosol list for dust1
          84             : integer :: idxdst2  = -1 ! index in aerosol list for dust2
          85             : integer :: idxdst3  = -1 ! index in aerosol list for dust3
          86             : integer :: idxdst4  = -1 ! index in aerosol list for dust4
          87             : integer :: idxbcphi = -1 ! index in aerosol list for Soot (BCPHIL)
          88             :
          89             : ! modal aerosols
          90             : logical :: clim_modal_aero = .false.
          91             : logical :: prog_modal_aero = .false.
          92             :
          93             : logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies
          94             :
          95             : integer, allocatable :: aer_cnst_idx(:,:)
          96             :
          97             : !===============================================================================
          98             : contains
          99             : !===============================================================================
         100             :
         101        1536 : subroutine nucleate_ice_cam_readnl(nlfile)
         102             :
         103             :   use namelist_utils,  only: find_group_name
         104             :   use units,           only: getunit, freeunit
         105             :   use spmd_utils,      only: mpicom, masterprocid, mpi_logical, mpi_real8
         106             :
         107             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         108             :
         109             :   ! Local variables
         110             :   integer :: unitn, ierr
         111             :   character(len=*), parameter :: subname = 'nucleate_ice_cam_readnl'
         112             :
         113             :   namelist /nucleate_ice_nl/ use_preexisting_ice, hist_preexisting_ice, &
         114             :        nucleate_ice_subgrid, nucleate_ice_subgrid_strat, nucleate_ice_strat, &
         115             :        nucleate_ice_incloud, nucleate_ice_use_troplev
         116             :
         117             :   !-----------------------------------------------------------------------------
         118             :
         119        1536 :   if (masterproc) then
         120           2 :      unitn = getunit()
         121           2 :      open( unitn, file=trim(nlfile), status='old' )
         122           2 :      call find_group_name(unitn, 'nucleate_ice_nl', status=ierr)
         123           2 :      if (ierr == 0) then
         124           2 :         read(unitn, nucleate_ice_nl, iostat=ierr)
         125           2 :         if (ierr /= 0) then
         126           0 :            call endrun(subname // ':: ERROR reading namelist')
         127             :         end if
         128             :      end if
         129           2 :      close(unitn)
         130           2 :      call freeunit(unitn)
         131             :
         132             :   end if
         133             :
         134             :   ! Broadcast namelist variables
         135        1536 :   call mpi_bcast(use_preexisting_ice,  1, mpi_logical,masterprocid, mpicom, ierr)
         136        1536 :   call mpi_bcast(hist_preexisting_ice, 1, mpi_logical,masterprocid, mpicom, ierr)
         137        1536 :   call mpi_bcast(nucleate_ice_subgrid, 1, mpi_real8,  masterprocid, mpicom, ierr)
         138        1536 :   call mpi_bcast(nucleate_ice_subgrid_strat, 1, mpi_real8,  masterprocid, mpicom, ierr)
         139        1536 :   call mpi_bcast(nucleate_ice_strat,   1, mpi_real8,  masterprocid, mpicom, ierr)
         140        1536 :   call mpi_bcast(nucleate_ice_incloud, 1, mpi_logical,masterprocid, mpicom, ierr)
         141        1536 :   call mpi_bcast(nucleate_ice_use_troplev, 1, mpi_logical,masterprocid, mpicom, ierr)
         142             :
         143        1536 : end subroutine nucleate_ice_cam_readnl
         144             :
         145             : !================================================================================================
         146             :
         147           0 : subroutine nucleate_ice_cam_register()
         148             :
         149             :    ! global scope for NAAI needed when clubb_do_icesuper=.true.
         150           0 :    call pbuf_add_field('NAAI',     'global', dtype_r8, (/pcols,pver/), naai_idx)
         151           0 :    call pbuf_add_field('NAAI_HOM', 'physpkg', dtype_r8, (/pcols,pver/), naai_hom_idx)
         152             :
         153           0 : end subroutine nucleate_ice_cam_register
         154             :
         155             : !================================================================================================
         156             :
         157           0 : subroutine nucleate_ice_cam_init(mincld_in, bulk_scale_in, pbuf2d, aero_props)
         158             :    use phys_control, only: phys_getopts
         159             :    use time_manager, only: is_first_step
         160             :
         161             :    real(r8), intent(in) :: mincld_in
         162             :    real(r8), intent(in) :: bulk_scale_in
         163             :    class(aerosol_properties), optional, intent(in) :: aero_props
         164             :
         165             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         166             :
         167             :    ! local variables
         168             :    integer :: iaer
         169             :    integer :: ierr
         170             :    integer :: ispc, ibin
         171             :    integer :: idxtmp
         172             :    integer :: nmodes
         173             :
         174             :    character(len=*), parameter :: routine = 'nucleate_ice_cam_init'
         175             :    logical :: history_cesm_forcing
         176             :
         177             :    character(len=32) :: tmpname
         178             :
         179             :    !--------------------------------------------------------------------------------------------
         180           0 :    call phys_getopts(prog_modal_aero_out = prog_modal_aero, history_cesm_forcing_out = history_cesm_forcing)
         181             :
         182           0 :    mincld     = mincld_in
         183           0 :    bulk_scale = bulk_scale_in
         184             :
         185           0 :    lq(:) = .false.
         186             :
         187           0 :    if (prog_modal_aero.and.use_preexisting_ice) then
         188             :
         189           0 :       if (.not. present(aero_props)) then
         190           0 :          call endrun(routine//' :  aero_props must be present')
         191             :       end if
         192             :
         193             :       ! constituent tendencies are calculated only if use_preexisting_ice is TRUE
         194             :       ! set lq for constituent tendencies --
         195             :
         196           0 :       allocate(aer_cnst_idx(aero_props%nbins(),0:maxval(aero_props%nspecies())), stat=ierr)
         197           0 :       if( ierr /= 0 ) then
         198           0 :          call endrun(routine//': aer_cnst_idx allocation failed')
         199             :       end if
         200           0 :       aer_cnst_idx = -1
         201             :
         202           0 :       do ibin = 1, aero_props%nbins()
         203           0 :          if (aero_props%icenuc_updates_num(ibin)) then
         204             :
         205             :             ! constituents of this bin will need to be updated
         206             :
         207           0 :             if (aero_props%icenuc_updates_mmr(ibin,0)) then ! species 0 indicates bin MMR
         208           0 :                call aero_props%amb_mmr_name( ibin, 0, tmpname)
         209             :             else
         210           0 :                call aero_props%amb_num_name( ibin, tmpname)
         211             :             end if
         212             :
         213           0 :             call cnst_get_ind(tmpname, idxtmp, abort=.false.)
         214           0 :             aer_cnst_idx(ibin,0) = idxtmp
         215           0 :             if (idxtmp>0) then
         216           0 :                lq(idxtmp) = .true.
         217             :             end if
         218             :
         219             :             ! iterate over the species within the bin
         220           0 :             do ispc = 1, aero_props%nspecies(ibin)
         221           0 :                if (aero_props%icenuc_updates_mmr(ibin,ispc)) then
         222             :                   ! this aerosol constituent will be updated
         223           0 :                   call aero_props%amb_mmr_name( ibin, ispc, tmpname)
         224           0 :                   call cnst_get_ind(tmpname, idxtmp, abort=.false.)
         225           0 :                   aer_cnst_idx(ibin,ispc) = idxtmp
         226           0 :                   if (idxtmp>0) then
         227           0 :                      lq(idxtmp) = .true.
         228             :                   end if
         229             :                end if
         230             :             end do
         231             :
         232             :          end if
         233             :       end do
         234             :
         235             :    end if
         236             :
         237             :    ! Initialize naai.
         238           0 :    if (is_first_step()) then
         239           0 :       call pbuf_set_field(pbuf2d, naai_idx, 0.0_r8)
         240             :    end if
         241             :
         242           0 :    if( masterproc ) then
         243           0 :       write(iulog,*) 'nucleate_ice parameters:'
         244           0 :       write(iulog,*) '  mincld                     = ', mincld_in
         245           0 :       write(iulog,*) '  bulk_scale                 = ', bulk_scale_in
         246           0 :       write(iulog,*) '  use_preexisiting_ice       = ', use_preexisting_ice
         247           0 :       write(iulog,*) '  hist_preexisiting_ice      = ', hist_preexisting_ice
         248           0 :       write(iulog,*) '  nucleate_ice_subgrid       = ', nucleate_ice_subgrid
         249           0 :       write(iulog,*) '  nucleate_ice_subgrid_strat = ', nucleate_ice_subgrid_strat
         250           0 :       write(iulog,*) '  nucleate_ice_strat         = ', nucleate_ice_strat
         251           0 :       write(iulog,*) '  nucleate_ice_incloud       = ', nucleate_ice_incloud
         252           0 :       write(iulog,*) '  nucleate_ice_use_troplev   = ', nucleate_ice_use_troplev
         253             :    end if
         254             :
         255           0 :    call cnst_get_ind('CLDLIQ', cldliq_idx)
         256           0 :    call cnst_get_ind('CLDICE', cldice_idx)
         257           0 :    call cnst_get_ind('NUMICE', numice_idx)
         258           0 :    qsatfac_idx  = pbuf_get_index('QSATFAC', ierr)
         259             :
         260           0 :    if (((nucleate_ice_subgrid .eq. -1._r8) .or. (nucleate_ice_subgrid_strat .eq. -1._r8)) .and. (qsatfac_idx .eq. -1)) then
         261           0 :      call endrun(routine//': ERROR qsatfac is required when subgrid = -1 or subgrid_strat = -1')
         262             :    end if
         263             :
         264           0 :    if (cam_physpkg_is("cam7")) then
         265             :       ! Updates for PUMAS v1.21+
         266           0 :       call addfld('NIHFTEN',  (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to homogenous freezing')
         267           0 :       call addfld('NIDEPTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to deposition nucleation')
         268           0 :       call addfld('NIIMMTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to immersion freezing')
         269           0 :       call addfld('NIMEYTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to meyers deposition')
         270             :    else
         271           0 :       call addfld('NIHF',  (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to homogenous freezing')
         272           0 :       call addfld('NIDEP', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to deposition nucleation')
         273           0 :       call addfld('NIIMM', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to immersion freezing')
         274           0 :       call addfld('NIMEY', (/ 'lev' /), 'A', '1/m3', 'Activated Ice Number Concentration due to meyers deposition')
         275             :    endif
         276             :
         277           0 :    call addfld('NIREGM',(/ 'lev' /), 'A', 'C', 'Ice Nucleation Temperature Threshold for Regime')
         278           0 :    call addfld('NISUBGRID',(/ 'lev' /), 'A', '', 'Ice Nucleation subgrid saturation factor')
         279           0 :    call addfld('NITROP_PD',(/ 'lev' /), 'A', '', 'Chemical Tropopause probability')
         280           0 :    if ( history_cesm_forcing ) then
         281           0 :       call add_default('NITROP_PD',8,' ')
         282             :    endif
         283             :
         284           0 :    if (use_preexisting_ice) then
         285           0 :       call addfld('fhom',      (/ 'lev' /), 'A','fraction', 'Fraction of cirrus where homogeneous freezing occur'   )
         286           0 :       call addfld ('WICE',     (/ 'lev' /), 'A','m/s','Vertical velocity Reduction caused by preexisting ice'  )
         287           0 :       call addfld ('WEFF',     (/ 'lev' /), 'A','m/s','Effective Vertical velocity for ice nucleation' )
         288             :
         289           0 :       if (cam_physpkg_is("cam7")) then
         290             :          ! Updates for PUMAS v1.21+
         291           0 :          call addfld ('INnso4TEN',   (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency so4 (in) to ice_nucleation')
         292           0 :          call addfld ('INnbcTEN',    (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency bc  (in) to ice_nucleation')
         293           0 :          call addfld ('INndustTEN',  (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency dust (in) ice_nucleation')
         294           0 :          call addfld ('INondustTEN',  (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency dust (out) from ice_nucleation')
         295             :          call addfld ('INhetTEN',    (/ 'lev' /), 'A','1/m3/s', &
         296           0 :               'Tendency for contribution for in-cloud ice number density increase by het nucleation in ice cloud')
         297             :          call addfld ('INhomTEN',    (/ 'lev' /), 'A','1/m3/s', &
         298           0 :               'Tendency for contribution for in-cloud ice number density increase by hom nucleation in ice cloud')
         299             :       else
         300           0 :          call addfld ('INnso4',   (/ 'lev' /), 'A','1/m3','Number Concentration so4 (in) to ice_nucleation')
         301           0 :          call addfld ('INnbc',    (/ 'lev' /), 'A','1/m3','Number Concentration bc  (in) to ice_nucleation')
         302           0 :          call addfld ('INndust',  (/ 'lev' /), 'A','1/m3','Number Concentration dust (in) ice_nucleation')
         303           0 :          call addfld ('INondust',  (/ 'lev' /), 'A','1/m3','Number Concentration dust (out) from ice_nucleation')
         304             :          call addfld ('INhet',    (/ 'lev' /), 'A','1/m3', &
         305           0 :               'contribution for in-cloud ice number density increase by het nucleation in ice cloud')
         306             :          call addfld ('INhom',    (/ 'lev' /), 'A','1/m3', &
         307           0 :               'contribution for in-cloud ice number density increase by hom nucleation in ice cloud')
         308             :       endif
         309             :
         310           0 :       call addfld ('INFrehom', (/ 'lev' /), 'A','frequency','hom IN frequency ice cloud')
         311           0 :       call addfld ('INFreIN',  (/ 'lev' /), 'A','frequency','frequency of ice nucleation occur')
         312             :
         313           0 :       if (hist_preexisting_ice) then
         314           0 :          call add_default ('WSUBI   ', 1, ' ')  ! addfld/outfld calls are in microp_aero
         315             :
         316           0 :          call add_default ('fhom    ', 1, ' ')
         317           0 :          call add_default ('WICE    ', 1, ' ')
         318           0 :          call add_default ('WEFF    ', 1, ' ')
         319           0 :          call add_default ('INnso4  ', 1, ' ')
         320           0 :          call add_default ('INnbc   ', 1, ' ')
         321           0 :          call add_default ('INndust ', 1, ' ')
         322           0 :          call add_default ('INhet   ', 1, ' ')
         323           0 :          call add_default ('INhom   ', 1, ' ')
         324           0 :          call add_default ('INFrehom', 1, ' ')
         325           0 :          call add_default ('INFreIN ', 1, ' ')
         326             :       end if
         327             :    end if
         328             :
         329             :    ! clim_modal_aero determines whether modal aerosols are used in the climate calculation.
         330             :    ! The modal aerosols can be either prognostic or prescribed.
         331           0 :    call rad_cnst_get_info(0, nmodes=nmodes)
         332             :
         333           0 :    clim_modal_aero = (nmodes > 0)
         334             :
         335           0 :    if (.not. clim_modal_aero) then
         336             :
         337             :       ! Props needed for BAM number concentration calcs.
         338             :
         339           0 :       call rad_cnst_get_info(0, naero=naer_all)
         340             :       allocate( &
         341           0 :          aername(naer_all),        &
         342           0 :          num_to_mass_aer(naer_all) )
         343             :
         344           0 :       do iaer = 1, naer_all
         345             :          call rad_cnst_get_aer_props(0, iaer, &
         346           0 :             aername         = aername(iaer), &
         347           0 :             num_to_mass_aer = num_to_mass_aer(iaer))
         348             :          ! Look for sulfate, dust, and soot in this list (Bulk aerosol only)
         349           0 :          if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
         350           0 :          if (trim(aername(iaer)) == 'DUST1') idxdst1 = iaer
         351           0 :          if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
         352           0 :          if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
         353           0 :          if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
         354           0 :          if (trim(aername(iaer)) == 'BCPHI') idxbcphi = iaer
         355             :       end do
         356             :    end if
         357             :
         358             :
         359             :    call nucleati_init(use_preexisting_ice, use_hetfrz_classnuc, nucleate_ice_incloud, iulog, pi, &
         360           0 :         mincld)
         361             :
         362             :    ! get indices for fields in the physics buffer
         363           0 :    ast_idx = pbuf_get_index('AST')
         364             :
         365           0 : end subroutine nucleate_ice_cam_init
         366             :
         367             : !================================================================================================
         368             :
         369           0 : subroutine nucleate_ice_cam_calc( &
         370           0 :    state, wsubi, pbuf, dtime, ptend, aero_props, aero_state )
         371             :
         372           0 :    use tropopause,     only: tropopause_findChemTrop
         373             :
         374             :    ! arguments
         375             :    type(physics_state), target, intent(in)    :: state
         376             :    real(r8),                    intent(in)    :: wsubi(:,:)
         377             :    type(physics_buffer_desc),   pointer       :: pbuf(:)
         378             :    real(r8),                    intent(in)    :: dtime
         379             :    type(physics_ptend),         intent(out)   :: ptend
         380             :    class(aerosol_properties),optional, intent(in) :: aero_props
         381             :    class(aerosol_state),optional, intent(in) :: aero_state
         382             :
         383             :    ! local workspace
         384             :
         385             :    ! naai and naai_hom are the outputs shared with the microphysics
         386           0 :    real(r8), pointer :: naai(:,:)       ! number of activated aerosol for ice nucleation
         387           0 :    real(r8), pointer :: naai_hom(:,:)   ! number of activated aerosol for ice nucleation (homogeneous freezing only)
         388             :
         389             :    integer :: lchnk, ncol
         390             :    integer :: itim_old
         391             :    integer :: i, k, l, m
         392             :
         393             :    character(len=32) :: spectype
         394             :
         395           0 :    real(r8), pointer :: t(:,:)          ! input temperature (K)
         396           0 :    real(r8), pointer :: qn(:,:)         ! input water vapor mixing ratio (kg/kg)
         397           0 :    real(r8), pointer :: qc(:,:)         ! cloud water mixing ratio (kg/kg)
         398           0 :    real(r8), pointer :: qi(:,:)         ! cloud ice mixing ratio (kg/kg)
         399           0 :    real(r8), pointer :: ni(:,:)         ! cloud ice number conc (1/kg)
         400           0 :    real(r8), pointer :: pmid(:,:)       ! pressure at layer midpoints (pa)
         401             :
         402           0 :    real(r8), pointer :: aer_mmr(:,:)    ! aerosol mass mixing ratio
         403             :
         404           0 :    real(r8), pointer :: ast(:,:)
         405             :    real(r8) :: icecldf(pcols,pver)  ! ice cloud fraction
         406           0 :    real(r8), pointer :: qsatfac(:,:)      ! Subgrid cloud water saturation scaling factor.
         407             :
         408             :    real(r8) :: rho(pcols,pver)      ! air density (kg m-3)
         409             :
         410           0 :    real(r8), allocatable :: naer2(:,:,:)    ! bulk aerosol number concentration (1/m3)
         411           0 :    real(r8), allocatable :: maerosol(:,:,:) ! bulk aerosol mass conc (kg/m3)
         412             :
         413             :    real(r8) :: qs(pcols)            ! liquid-ice weighted sat mixing rat (kg/kg)
         414             :    real(r8) :: es(pcols)            ! liquid-ice weighted sat vapor press (pa)
         415             :    real(r8) :: gammas(pcols)        ! parameter for cond/evap of cloud water
         416             :    integer  :: troplev(pcols)       ! tropopause level
         417             :
         418             :    real(r8) :: relhum(pcols,pver)  ! relative humidity
         419             :    real(r8) :: icldm(pcols,pver)   ! ice cloud fraction
         420             :
         421             :    real(r8) :: dst_num                               ! total dust aerosol number (#/cm^3)
         422             :    real(r8) :: dso4_num                               ! so4 aerosol number (#/cm^3)
         423             :    real(r8) :: so4_num                               ! so4 aerosol number (#/cm^3)
         424             :    real(r8) :: soot_num                              ! soot (hydrophilic) aerosol number (#/cm^3)
         425             :    real(r8) :: wght
         426             :    real(r8) :: oso4_num
         427             :    real(r8) :: odst_num
         428             :    real(r8) :: osoot_num
         429             :    real(r8) :: so4_num_st_cr_tot
         430             :    real(r8) :: ramp
         431             :
         432             :    real(r8) :: subgrid(pcols,pver)
         433             :    real(r8) :: trop_pd(pcols,pver)
         434             :
         435             :    ! For pre-existing ice
         436             :    real(r8) :: fhom(pcols,pver)    ! how much fraction of cloud can reach Shom
         437             :    real(r8) :: wice(pcols,pver)    ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at Shom
         438             :    real(r8) :: weff(pcols,pver)    ! effective Vertical velocity for ice nucleation (m/s); weff=wsubi-wice
         439             :    real(r8) :: INnso4(pcols,pver)   ! #/m3, so4 aerosol number used for ice nucleation
         440             :    real(r8) :: INnbc(pcols,pver)    ! #/m3, bc aerosol number used for ice nucleation
         441             :    real(r8) :: INndust(pcols,pver)  ! #/m3, dust aerosol number used for ice nucleation
         442             :    real(r8) :: INondust(pcols,pver)  ! #/m3, dust aerosol number used for ice nucleation
         443             :    real(r8) :: INhet(pcols,pver)    ! #/m3, ice number from het freezing
         444             :    real(r8) :: INhom(pcols,pver)    ! #/m3, ice number from hom freezing
         445             :    real(r8) :: INFrehom(pcols,pver) !  hom freezing occurence frequency.  1 occur, 0 not occur.
         446             :    real(r8) :: INFreIN(pcols,pver)  !  ice nucleation occerence frequency.   1 occur, 0 not occur.
         447             :
         448             :    ! history output for ice nucleation
         449             :    real(r8) :: nihf(pcols,pver)  !output number conc of ice nuclei due to heterogenous freezing (1/m3)
         450             :    real(r8) :: niimm(pcols,pver) !output number conc of ice nuclei due to immersion freezing (hetero nuc) (1/m3)
         451             :    real(r8) :: nidep(pcols,pver) !output number conc of ice nuclei due to deoposion nucleation (hetero nuc) (1/m3)
         452             :    real(r8) :: nimey(pcols,pver) !output number conc of ice nuclei due to meyers deposition (1/m3)
         453             :    real(r8) :: regm(pcols,pver)  !output temperature thershold for nucleation regime
         454             :
         455             :    real(r8) :: size_wghts(pcols,pver)
         456             :    real(r8) :: type_wghts(pcols,pver)
         457             :    real(r8), pointer :: num_col(:,:)
         458             :    real(r8) :: dust_num_col(pcols,pver)
         459             :    real(r8) :: sulf_num_col(pcols,pver)
         460             :    real(r8) :: soot_num_col(pcols,pver)
         461             :    real(r8) :: sulf_num_tot_col(pcols,pver)
         462             :
         463             :    integer :: idxtmp
         464           0 :    real(r8), pointer :: amb_num(:,:)
         465           0 :    real(r8), pointer :: amb_mmr(:,:)
         466           0 :    real(r8), pointer :: cld_num(:,:)
         467           0 :    real(r8), pointer :: cld_mmr(:,:)
         468             :
         469             :    real(r8) :: delmmr, delmmr_sum
         470             :    real(r8) :: delnum, delnum_sum
         471             :
         472             :    real(r8), parameter :: per_cm3 = 1.e-6_r8 ! factor for m-3 to cm-3 conversions
         473             :
         474             :    !-------------------------------------------------------------------------------
         475             :
         476           0 :    lchnk = state%lchnk
         477           0 :    ncol  = state%ncol
         478           0 :    t     => state%t
         479           0 :    qn    => state%q(:,:,1)
         480           0 :    qc    => state%q(:,:,cldliq_idx)
         481           0 :    qi    => state%q(:,:,cldice_idx)
         482           0 :    ni    => state%q(:,:,numice_idx)
         483           0 :    pmid  => state%pmid
         484             :
         485           0 :    rho(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:))
         486             :
         487           0 :    if (clim_modal_aero) then
         488             :
         489           0 :       call physics_ptend_init(ptend, state%psetcols, 'nucleatei', lq=lq)
         490             :
         491             :    else
         492             :       ! init number/mass arrays for bulk aerosols
         493             :       allocate( &
         494             :            naer2(pcols,pver,naer_all), &
         495           0 :            maerosol(pcols,pver,naer_all))
         496             :
         497           0 :       do m = 1, naer_all
         498           0 :          call rad_cnst_get_aer_mmr(0, m, state, pbuf, aer_mmr)
         499           0 :          maerosol(:ncol,:,m) = aer_mmr(:ncol,:)*rho(:ncol,:)
         500             :
         501           0 :          if (m .eq. idxsul) then
         502           0 :             naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)*bulk_scale
         503             :          else
         504           0 :             naer2(:ncol,:,m) = maerosol(:ncol,:,m)*num_to_mass_aer(m)
         505             :          end if
         506             :       end do
         507             :
         508           0 :       call physics_ptend_init(ptend, state%psetcols, 'nucleatei')
         509             :    end if
         510             :
         511           0 :    itim_old = pbuf_old_tim_idx()
         512           0 :    call pbuf_get_field(pbuf, ast_idx, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         513             :
         514           0 :    icecldf(:ncol,:pver) = ast(:ncol,:pver)
         515             :
         516             :    ! naai and naai_hom are the outputs from this parameterization
         517           0 :    call pbuf_get_field(pbuf, naai_idx, naai)
         518           0 :    call pbuf_get_field(pbuf, naai_hom_idx, naai_hom)
         519           0 :    naai(1:ncol,1:pver)     = 0._r8
         520           0 :    naai_hom(1:ncol,1:pver) = 0._r8
         521             :
         522             :    ! Use the same criteria that is used in chemistry and in CLUBB (for cloud fraction)
         523             :    ! to determine whether to use tropospheric or stratospheric settings. Include the
         524             :    ! tropopause level so that the cold point tropopause will use the stratospheric values.
         525             :    !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         526           0 :    troplev(:) = 0
         527             :    !REMOVECAM_END
         528           0 :    call tropopause_findChemTrop(state, troplev)
         529             :
         530           0 :    if ((nucleate_ice_subgrid .eq. -1._r8) .or. (nucleate_ice_subgrid_strat .eq. -1._r8)) then
         531           0 :       call pbuf_get_field(pbuf, qsatfac_idx, qsatfac)
         532             :    end if
         533             :
         534           0 :    trop_pd(:,:) = 0._r8
         535             :
         536           0 :    do k = top_lev, pver
         537           0 :       do i = 1, ncol
         538           0 :          trop_pd(i, troplev(i)) = 1._r8
         539             :
         540           0 :          if (k <= troplev(i)) then
         541           0 :             if (nucleate_ice_subgrid_strat .eq. -1._r8) then
         542           0 :                subgrid(i, k) = 1._r8 / qsatfac(i, k)
         543             :             else
         544           0 :                subgrid(i, k) = nucleate_ice_subgrid_strat
         545             :             end if
         546             :          else
         547           0 :             if (nucleate_ice_subgrid .eq. -1._r8) then
         548           0 :                subgrid(i, k) = 1._r8 / qsatfac(i, k)
         549             :             else
         550           0 :                subgrid(i, k) = nucleate_ice_subgrid
         551             :             end if
         552             :          end if
         553             :       end do
         554             :    end do
         555             :
         556             :
         557             :    ! initialize history output fields for ice nucleation
         558           0 :    nihf(1:ncol,1:pver)  = 0._r8
         559           0 :    niimm(1:ncol,1:pver) = 0._r8
         560           0 :    nidep(1:ncol,1:pver) = 0._r8
         561           0 :    nimey(1:ncol,1:pver) = 0._r8
         562             :
         563           0 :    regm(1:ncol,1:pver) = 0._r8
         564             :
         565           0 :    if (use_preexisting_ice) then
         566           0 :       fhom(:,:)     = 0.0_r8
         567           0 :       wice(:,:)     = 0.0_r8
         568           0 :       weff(:,:)     = 0.0_r8
         569           0 :       INnso4(:,:)   = 0.0_r8
         570           0 :       INnbc(:,:)    = 0.0_r8
         571           0 :       INndust(:,:)  = 0.0_r8
         572           0 :       INondust(:,:)  = 0.0_r8
         573           0 :       INhet(:,:)    = 0.0_r8
         574           0 :       INhom(:,:)    = 0.0_r8
         575           0 :       INFrehom(:,:) = 0.0_r8
         576           0 :       INFreIN(:,:)  = 0.0_r8
         577             :    endif
         578             :
         579           0 :    do k = top_lev, pver
         580             :       ! Get humidity and saturation vapor pressures
         581           0 :       call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gammas(1:ncol))
         582             :
         583           0 :       do i = 1, ncol
         584             :
         585           0 :          relhum(i,k) = qn(i,k)/qs(i)
         586             :
         587             :          ! get cloud fraction, check for minimum
         588           0 :          icldm(i,k) = max(icecldf(i,k), mincld)
         589             :
         590             :       end do
         591             :    end do
         592             :
         593           0 :    dust_num_col = 0._r8
         594           0 :    sulf_num_col = 0._r8
         595           0 :    sulf_num_tot_col = 0._r8
         596           0 :    soot_num_col = 0._r8
         597             :
         598           0 :    if (clim_modal_aero) then
         599             :
         600           0 :       if (.not.(present(aero_props).and.present(aero_state))) then
         601           0 :          call endrun('nucleate_ice_cam_calc: aero_props and aero_state must be present')
         602             :       end if
         603             :
         604             :       ! collect number densities (#/cm^3) for dust, sulfate, and soot
         605             :       call aero_state%nuclice_get_numdens( aero_props, use_preexisting_ice, ncol, pver, rho, &
         606           0 :                                            dust_num_col, sulf_num_col, soot_num_col, sulf_num_tot_col )
         607             :
         608             :    else
         609             :       ! for bulk model
         610           0 :       dust_num_col(:ncol,:) = naer2(:ncol,:,idxdst1)/25._r8 * per_cm3 & ! #/cm3
         611           0 :                             + naer2(:ncol,:,idxdst2)/25._r8 * per_cm3 &
         612           0 :                             + naer2(:ncol,:,idxdst3)/25._r8 * per_cm3 &
         613           0 :                             + naer2(:ncol,:,idxdst4)/25._r8 * per_cm3
         614           0 :       sulf_num_col(:ncol,:) = naer2(:ncol,:,idxsul)/25._r8 * per_cm3
         615           0 :       soot_num_col(:ncol,:) = naer2(:ncol,:,idxbcphi)/25._r8 * per_cm3
         616             :    endif
         617             :
         618           0 :    kloop: do k = top_lev, pver
         619           0 :       iloop: do i = 1, ncol
         620             :
         621           0 :          so4_num_st_cr_tot = 0._r8
         622             :
         623           0 :          freezing: if (t(i,k) < tmelt - 5._r8) then
         624             :
         625             :             ! set aerosol number for so4, soot, and dust with units #/cm^3
         626           0 :             so4_num = sulf_num_col(i,k)
         627           0 :             dst_num = dust_num_col(i,k)
         628           0 :             so4_num_st_cr_tot=sulf_num_tot_col(i,k)
         629             :
         630             :             ! *** Turn off soot nucleation ***
         631           0 :             soot_num = 0.0_r8
         632             :
         633           0 :             if (cam_physpkg_is("cam7")) then
         634             :
         635             :                call nucleati( &
         636           0 :                     wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k),   &
         637           0 :                     qc(i,k), qi(i,k), ni(i,k), rho(i,k),                      &
         638             :                     so4_num, dst_num, soot_num, subgrid(i,k),                 &
         639           0 :                     naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
         640             :                     wice(i,k), weff(i,k), fhom(i,k), regm(i,k),               &
         641             :                     oso4_num, odst_num, osoot_num, &
         642           0 :                     call_frm_zm_in = .false., add_preexisting_ice_in = .false.)
         643             :
         644             :             else
         645             :
         646             :                call nucleati( &
         647           0 :                     wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k),   &
         648           0 :                     qc(i,k), qi(i,k), ni(i,k), rho(i,k),                      &
         649             :                     so4_num, dst_num, soot_num, subgrid(i,k),                 &
         650           0 :                     naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), &
         651             :                     wice(i,k), weff(i,k), fhom(i,k), regm(i,k),               &
         652           0 :                     oso4_num, odst_num, osoot_num)
         653             :
         654             :             end if
         655             :
         656             :             ! Move aerosol used for nucleation from interstial to cloudborne,
         657             :             ! otherwise the same coarse mode aerosols will be available again
         658             :             ! in the next timestep and will supress homogeneous freezing.
         659             :
         660             :
         661           0 :             if (prog_modal_aero .and. use_preexisting_ice) then
         662             :
         663             :                ! compute tendencies for transported aerosol constituents
         664             :                ! and update not-transported constituents
         665             :
         666           0 :                do m = 1, aero_props%nbins()
         667             :
         668           0 :                   if (aero_props%icenuc_updates_num(m)) then
         669             :
         670             :                      ! constituents of this bin will need to be updated
         671             :
         672           0 :                      call aero_state%get_ambient_num(m, amb_num)
         673           0 :                      call aero_state%get_cldbrne_num(m, cld_num)
         674             :
         675           0 :                      if (amb_num(i,k)>0._r8) then
         676           0 :                         delmmr_sum = 0._r8
         677           0 :                         delnum_sum = 0._r8
         678             :
         679             :                         ! iterate over the species within the bin
         680           0 :                         do l = 1, aero_props%nspecies(m)
         681           0 :                            if (aero_props%icenuc_updates_mmr(m,l)) then
         682             :
         683           0 :                               call aero_props%species_type(m, l, spectype)
         684           0 :                               call aero_state%icenuc_size_wght( m, i,k, spectype, use_preexisting_ice, wght)
         685             :
         686           0 :                               if (wght>0._r8) then
         687             :
         688             :                                  ! this aerosol constituent will be updated
         689             :
         690           0 :                                  idxtmp = aer_cnst_idx(m,l)
         691             :
         692           0 :                                  call aero_state%get_ambient_mmr(l,m,amb_mmr)
         693           0 :                                  call aero_state%get_cldbrne_mmr(l,m,cld_mmr)
         694             :
         695             :                                  ! determine change in aerosol mass
         696           0 :                                  delmmr = 0._r8
         697           0 :                                  delnum = 0._r8
         698           0 :                                  if (trim(spectype)=='dust') then
         699           0 :                                     if (dst_num>0._r8) then
         700           0 :                                        delmmr = (odst_num / dst_num) * icldm(i,k) * amb_mmr(i,k) * wght
         701           0 :                                        delnum = (odst_num * icldm(i,k)) /rho(i,k)/per_cm3
         702             :                                     endif
         703           0 :                                  elseif (trim(spectype)=='sulfate') then
         704           0 :                                     if (so4_num>0._r8) then
         705           0 :                                        delmmr = (oso4_num / so4_num) * icldm(i,k) * amb_mmr(i,k) * wght
         706           0 :                                        delnum = (oso4_num * icldm(i,k)) /rho(i,k)/per_cm3
         707             :                                     endif
         708             :                                  endif
         709             :
         710           0 :                                  if (idxtmp>0) then
         711             :                                     ! constituent tendency (for transported species)
         712           0 :                                     ptend%q(i,k,idxtmp) = -delmmr/dtime
         713             :                                  else
         714             :                                     ! apply change of mass to not-transported species
         715           0 :                                     amb_mmr(i,k) = amb_mmr(i,k) - delmmr
         716             :                                  endif
         717           0 :                                  cld_mmr(i,k) = cld_mmr(i,k) + delmmr
         718             :
         719           0 :                                  delmmr_sum = delmmr_sum + delmmr
         720           0 :                                  delnum_sum = delnum_sum + delnum
         721             :                               end if
         722             :                            end if
         723             :                         end do
         724             :
         725           0 :                         idxtmp = aer_cnst_idx(m,0)
         726             :
         727             :                         ! update aerosol state bin and tendency for grid box i,k
         728           0 :                         call aero_state%update_bin( m,i,k, delmmr_sum, delnum_sum, idxtmp, dtime, ptend%q )
         729             :
         730             :                      end if
         731             :
         732             :                   end if
         733             :                end do
         734             :
         735             :             end if
         736             :
         737             :
         738             :             ! Liu&Penner does not generate enough nucleation in the polar winter
         739             :             ! stratosphere, which affects surface area density, dehydration and
         740             :             ! ozone chemistry. Part of this is that there are a larger number of
         741             :             ! particles in the accumulation mode than in the Aitken mode. In volcanic
         742             :             ! periods, the coarse mode may also be important. As a short
         743             :             ! term work around, include the accumulation and coarse mode particles
         744             :             ! and assume a larger fraction of the sulfates nucleate in the polar
         745             :             ! stratosphere.
         746             :             !
         747             :             ! Do not include the tropopause level, as stratospheric aerosols
         748             :             ! only exist above the tropopause level.
         749             :             !
         750             :             ! NOTE: This may still not represent the proper particles that
         751             :             ! participate in nucleation, because it doesn't include STS and NAT
         752             :             ! particles. It may not represent the proper saturation threshold for
         753             :             ! nucleation, and wsubi from CLUBB is probably not representative of
         754             :             ! wave driven varaibility in the polar stratosphere.
         755           0 :             if (nucleate_ice_use_troplev .and. clim_modal_aero) then
         756           0 :                if ((k < troplev(i)) .and. (nucleate_ice_strat > 0._r8) .and. (oso4_num > 0._r8)) then
         757           0 :                   dso4_num = max(0._r8, (nucleate_ice_strat*so4_num_st_cr_tot - oso4_num) * 1e6_r8 / rho(i,k))
         758           0 :                   naai(i,k) = naai(i,k) + dso4_num
         759           0 :                   nihf(i,k) = nihf(i,k) + dso4_num
         760             :                endif
         761             :             else
         762             :                ! This maintains backwards compatibility with the previous version.
         763           0 :                if (pmid(i,k) <= 12500._r8 .and. pmid(i,k) > 100._r8 .and. abs(state%lat(i)) >= 60._r8 * pi / 180._r8) then
         764           0 :                   ramp = 1._r8 - min(1._r8, max(0._r8, (pmid(i,k) - 10000._r8) / 2500._r8))
         765             :
         766           0 :                   if (oso4_num > 0._r8) then
         767           0 :                      dso4_num = (max(oso4_num, ramp * nucleate_ice_strat * so4_num) - oso4_num) * 1e6_r8 / rho(i,k)
         768           0 :                      naai(i,k) = naai(i,k) + dso4_num
         769           0 :                      nihf(i,k) = nihf(i,k) + dso4_num
         770             :                   end if
         771             :                end if
         772             :             end if
         773             :
         774           0 :             if (cam_physpkg_is("cam7")) then
         775             :                !Updates for pumas v1.21+
         776             :
         777           0 :                naai_hom(i,k) = nihf(i,k)/dtime
         778           0 :                naai(i,k)= naai(i,k)/dtime
         779             :
         780             :                ! output activated ice (convert from #/kg -> #/m3/s)
         781           0 :                nihf(i,k)     = nihf(i,k) *rho(i,k)/dtime
         782           0 :                niimm(i,k)    = niimm(i,k)*rho(i,k)/dtime
         783           0 :                nidep(i,k)    = nidep(i,k)*rho(i,k)/dtime
         784           0 :                nimey(i,k)    = nimey(i,k)*rho(i,k)/dtime
         785             :
         786           0 :                if (use_preexisting_ice) then
         787           0 :                   INnso4(i,k) =so4_num*1e6_r8/dtime  ! (convert from #/cm3 -> #/m3/s)
         788           0 :                   INnbc(i,k)  =soot_num*1e6_r8/dtime
         789           0 :                   INndust(i,k)=dst_num*1e6_r8/dtime
         790           0 :                   INondust(i,k)=odst_num*1e6_r8/dtime
         791           0 :                   INFreIN(i,k)=1.0_r8          ! 1,ice nucleation occur
         792           0 :                   INhet(i,k) = (niimm(i,k) + nidep(i,k))   ! #/m3/s, nimey not in cirrus
         793           0 :                   INhom(i,k) = nihf(i,k)                 ! #/m3/s
         794           0 :                   if (INhom(i,k).gt.1e3_r8)   then ! > 1/L
         795           0 :                      INFrehom(i,k)=1.0_r8       ! 1, hom freezing occur
         796             :                   endif
         797             :
         798             :                   ! exclude  no ice nucleaton
         799           0 :                   if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8))   then
         800           0 :                      INnso4(i,k) =0.0_r8
         801           0 :                      INnbc(i,k)  =0.0_r8
         802           0 :                      INndust(i,k)=0.0_r8
         803           0 :                      INondust(i,k)=0.0_r8
         804           0 :                      INFreIN(i,k)=0.0_r8
         805           0 :                      INhet(i,k) = 0.0_r8
         806           0 :                      INhom(i,k) = 0.0_r8
         807           0 :                      INFrehom(i,k)=0.0_r8
         808           0 :                      wice(i,k) = 0.0_r8
         809           0 :                      weff(i,k) = 0.0_r8
         810           0 :                      fhom(i,k) = 0.0_r8
         811             :                   endif
         812             :                endif
         813             :
         814             :             else ! Not cam7
         815             :
         816           0 :                naai_hom(i,k) = nihf(i,k)
         817             :
         818             :                ! output activated ice (convert from #/kg -> #/m3/s)
         819           0 :                nihf(i,k)     = nihf(i,k) *rho(i,k)
         820           0 :                niimm(i,k)    = niimm(i,k)*rho(i,k)
         821           0 :                nidep(i,k)    = nidep(i,k)*rho(i,k)
         822           0 :                nimey(i,k)    = nimey(i,k)*rho(i,k)
         823             :
         824           0 :                if (use_preexisting_ice) then
         825           0 :                   INnso4(i,k) =so4_num*1e6_r8 ! (convert from #/cm3 -> #/m3/s)
         826           0 :                   INnbc(i,k)  =soot_num*1e6_r8
         827           0 :                   INndust(i,k)=dst_num*1e6_r8
         828           0 :                   INondust(i,k)=odst_num*1e6_r8
         829           0 :                   INFreIN(i,k)=1.0_r8          ! 1,ice nucleation occur
         830           0 :                   INhet(i,k) = (niimm(i,k) + nidep(i,k))   ! #/m3, nimey not in cirrus
         831           0 :                   INhom(i,k) = nihf(i,k)                 ! #/m3
         832           0 :                   if (INhom(i,k).gt.1e3_r8)   then ! > 1/L
         833           0 :                      INFrehom(i,k)=1.0_r8       ! 1, hom freezing occur
         834             :                   endif
         835             :
         836             :                   ! exclude  no ice nucleaton
         837           0 :                   if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8))   then
         838           0 :                      INnso4(i,k) =0.0_r8
         839           0 :                      INnbc(i,k)  =0.0_r8
         840           0 :                      INndust(i,k)=0.0_r8
         841           0 :                      INondust(i,k)=0.0_r8
         842           0 :                      INFreIN(i,k)=0.0_r8
         843           0 :                      INhet(i,k) = 0.0_r8
         844           0 :                      INhom(i,k) = 0.0_r8
         845           0 :                      INFrehom(i,k)=0.0_r8
         846           0 :                      wice(i,k) = 0.0_r8
         847           0 :                      weff(i,k) = 0.0_r8
         848           0 :                      fhom(i,k) = 0.0_r8
         849             :                   endif
         850             :                end if
         851             :
         852             :             end if ! cam7
         853             :          end if freezing
         854             :       end do iloop
         855             :    end do kloop
         856             :
         857           0 :    if (.not. clim_modal_aero) then
         858           0 :       deallocate( &
         859             :            naer2, &
         860           0 :            maerosol)
         861             :    end if
         862             :
         863           0 :    if (cam_physpkg_is("cam7")) then
         864             :       ! Updates for PUMAS v1.21+
         865           0 :       call outfld('NIHFTEN',   nihf, pcols, lchnk)
         866           0 :       call outfld('NIIMMTEN', niimm, pcols, lchnk)
         867           0 :       call outfld('NIDEPTEN', nidep, pcols, lchnk)
         868           0 :       call outfld('NIMEYTEN', nimey, pcols, lchnk)
         869             :    else
         870           0 :       call outfld('NIHF',   nihf, pcols, lchnk)
         871           0 :       call outfld('NIIMM', niimm, pcols, lchnk)
         872           0 :       call outfld('NIDEP', nidep, pcols, lchnk)
         873           0 :       call outfld('NIMEY', nimey, pcols, lchnk)
         874             :    end if
         875           0 :    call outfld('NIREGM', regm, pcols, lchnk)
         876           0 :    call outfld('NISUBGRID', subgrid, pcols, lchnk)
         877           0 :    call outfld('NITROP_PD', trop_pd, pcols, lchnk)
         878             :
         879           0 :    if (use_preexisting_ice) then
         880           0 :       call outfld( 'fhom' , fhom, pcols, lchnk)
         881           0 :       call outfld( 'WICE' , wice, pcols, lchnk)
         882           0 :       call outfld( 'WEFF' , weff, pcols, lchnk)
         883           0 :       if (cam_physpkg_is("cam7")) then
         884             :          ! Updates for PUMAS v1.21+
         885           0 :          call outfld('INnso4TEN',INnso4 , pcols,lchnk)
         886           0 :          call outfld('INnbcTEN',INnbc  , pcols,lchnk)
         887           0 :          call outfld('INndustTEN',INndust, pcols,lchnk)
         888           0 :          call outfld('INondustTEN',INondust, pcols,lchnk)
         889           0 :          call outfld('INhetTEN',INhet  , pcols,lchnk)
         890           0 :          call outfld('INhomTEN',INhom  , pcols,lchnk)
         891             :       else
         892           0 :          call outfld('INnso4  ',INnso4 , pcols,lchnk)
         893           0 :          call outfld('INnbc   ',INnbc  , pcols,lchnk)
         894           0 :          call outfld('INndust ',INndust, pcols,lchnk)
         895           0 :          call outfld('INondust ',INondust, pcols,lchnk)
         896           0 :          call outfld('INhet   ',INhet  , pcols,lchnk)
         897           0 :          call outfld('INhom   ',INhom  , pcols,lchnk)
         898             :       end if
         899           0 :       call outfld('INFrehom',INFrehom,pcols,lchnk)
         900           0 :       call outfld('INFreIN ',INFreIN, pcols,lchnk)
         901             :    end if
         902             :
         903           0 : end subroutine nucleate_ice_cam_calc
         904             :
         905             : !================================================================================================
         906             :
         907             : end module nucleate_ice_cam

Generated by: LCOV version 1.14