LCOV - code coverage report
Current view: | top level - physics/cam - nucleate_ice_cam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 271 | 383 | 70.8 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 4 | 100.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     4469064 : 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        1536 : subroutine nucleate_ice_cam_register()
         148             :
         149             :    ! global scope for NAAI needed when clubb_do_icesuper=.true.
         150        1536 :    call pbuf_add_field('NAAI',     'global', dtype_r8, (/pcols,pver/), naai_idx)
         151        1536 :    call pbuf_add_field('NAAI_HOM', 'physpkg', dtype_r8, (/pcols,pver/), naai_hom_idx)
         152             :
         153        1536 : end subroutine nucleate_ice_cam_register
         154             :
         155             : !================================================================================================
         156             :
         157        1536 : 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        1536 :    call phys_getopts(prog_modal_aero_out = prog_modal_aero, history_cesm_forcing_out = history_cesm_forcing)
         181             :
         182        1536 :    mincld     = mincld_in
         183        1536 :    bulk_scale = bulk_scale_in
         184             :
         185        1536 :    lq(:) = .false.
         186             :
         187        1536 :    if (prog_modal_aero.and.use_preexisting_ice) then
         188             :
         189        1536 :       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       13824 :       allocate(aer_cnst_idx(aero_props%nbins(),0:maxval(aero_props%nspecies())), stat=ierr)
         197        1536 :       if( ierr /= 0 ) then
         198           0 :          call endrun(routine//': aer_cnst_idx allocation failed')
         199             :       end if
         200       55296 :       aer_cnst_idx = -1
         201             :
         202        9216 :       do ibin = 1, aero_props%nbins()
         203        7680 :          if (aero_props%icenuc_updates_num(ibin)) then
         204             :
         205             :             ! constituents of this bin will need to be updated
         206             :
         207        1536 :             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        1536 :                call aero_props%amb_num_name( ibin, tmpname)
         211             :             end if
         212             :
         213        1536 :             call cnst_get_ind(tmpname, idxtmp, abort=.false.)
         214        1536 :             aer_cnst_idx(ibin,0) = idxtmp
         215        1536 :             if (idxtmp>0) then
         216        1536 :                lq(idxtmp) = .true.
         217             :             end if
         218             :
         219             :             ! iterate over the species within the bin
         220        6144 :             do ispc = 1, aero_props%nspecies(ibin)
         221        6144 :                if (aero_props%icenuc_updates_mmr(ibin,ispc)) then
         222             :                   ! this aerosol constituent will be updated
         223        1536 :                   call aero_props%amb_mmr_name( ibin, ispc, tmpname)
         224        1536 :                   call cnst_get_ind(tmpname, idxtmp, abort=.false.)
         225        1536 :                   aer_cnst_idx(ibin,ispc) = idxtmp
         226        1536 :                   if (idxtmp>0) then
         227        1536 :                      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        1536 :    if (is_first_step()) then
         239         768 :       call pbuf_set_field(pbuf2d, naai_idx, 0.0_r8)
         240             :    end if
         241             :
         242        1536 :    if( masterproc ) then
         243           2 :       write(iulog,*) 'nucleate_ice parameters:'
         244           2 :       write(iulog,*) '  mincld                     = ', mincld_in
         245           2 :       write(iulog,*) '  bulk_scale                 = ', bulk_scale_in
         246           2 :       write(iulog,*) '  use_preexisiting_ice       = ', use_preexisting_ice
         247           2 :       write(iulog,*) '  hist_preexisiting_ice      = ', hist_preexisting_ice
         248           2 :       write(iulog,*) '  nucleate_ice_subgrid       = ', nucleate_ice_subgrid
         249           2 :       write(iulog,*) '  nucleate_ice_subgrid_strat = ', nucleate_ice_subgrid_strat
         250           2 :       write(iulog,*) '  nucleate_ice_strat         = ', nucleate_ice_strat
         251           2 :       write(iulog,*) '  nucleate_ice_incloud       = ', nucleate_ice_incloud
         252           2 :       write(iulog,*) '  nucleate_ice_use_troplev   = ', nucleate_ice_use_troplev
         253             :    end if
         254             :
         255        1536 :    call cnst_get_ind('CLDLIQ', cldliq_idx)
         256        1536 :    call cnst_get_ind('CLDICE', cldice_idx)
         257        1536 :    call cnst_get_ind('NUMICE', numice_idx)
         258        1536 :    qsatfac_idx  = pbuf_get_index('QSATFAC', ierr)
         259             :
         260        1536 :    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        1536 :    if (cam_physpkg_is("cam7")) then
         265             :       ! Updates for PUMAS v1.21+
         266        3072 :       call addfld('NIHFTEN',  (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to homogenous freezing')
         267        3072 :       call addfld('NIDEPTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to deposition nucleation')
         268        3072 :       call addfld('NIIMMTEN', (/ 'lev' /), 'A', '1/m3/s', 'Activated Ice Number Concentration tendency due to immersion freezing')
         269        3072 :       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        3072 :    call addfld('NIREGM',(/ 'lev' /), 'A', 'C', 'Ice Nucleation Temperature Threshold for Regime')
         278        3072 :    call addfld('NISUBGRID',(/ 'lev' /), 'A', '', 'Ice Nucleation subgrid saturation factor')
         279        3072 :    call addfld('NITROP_PD',(/ 'lev' /), 'A', '', 'Chemical Tropopause probability')
         280        1536 :    if ( history_cesm_forcing ) then
         281           0 :       call add_default('NITROP_PD',8,' ')
         282             :    endif
         283             :
         284        1536 :    if (use_preexisting_ice) then
         285        3072 :       call addfld('fhom',      (/ 'lev' /), 'A','fraction', 'Fraction of cirrus where homogeneous freezing occur'   )
         286        3072 :       call addfld ('WICE',     (/ 'lev' /), 'A','m/s','Vertical velocity Reduction caused by preexisting ice'  )
         287        3072 :       call addfld ('WEFF',     (/ 'lev' /), 'A','m/s','Effective Vertical velocity for ice nucleation' )
         288             :
         289        1536 :       if (cam_physpkg_is("cam7")) then
         290             :          ! Updates for PUMAS v1.21+
         291        3072 :          call addfld ('INnso4TEN',   (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency so4 (in) to ice_nucleation')
         292        3072 :          call addfld ('INnbcTEN',    (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency bc  (in) to ice_nucleation')
         293        3072 :          call addfld ('INndustTEN',  (/ 'lev' /), 'A','1/m3/s','Number Concentration tendency dust (in) ice_nucleation')
         294        3072 :          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        3072 :               '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        3072 :               '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        3072 :       call addfld ('INFrehom', (/ 'lev' /), 'A','frequency','hom IN frequency ice cloud')
         311        3072 :       call addfld ('INFreIN',  (/ 'lev' /), 'A','frequency','frequency of ice nucleation occur')
         312             :
         313        1536 :       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        1536 :    call rad_cnst_get_info(0, nmodes=nmodes)
         332             :
         333        1536 :    clim_modal_aero = (nmodes > 0)
         334             :
         335        1536 :    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        1536 :         mincld)
         361             :
         362             :    ! get indices for fields in the physics buffer
         363        1536 :    ast_idx = pbuf_get_index('AST')
         364             :
         365        1536 : end subroutine nucleate_ice_cam_init
         366             :
         367             : !================================================================================================
         368             :
         369   192103704 : subroutine nucleate_ice_cam_calc( &
         370     4467528 :    state, wsubi, pbuf, dtime, ptend, aero_props, aero_state )
         371             :
         372        1536 :    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     4467528 :    real(r8), pointer :: naai(:,:)       ! number of activated aerosol for ice nucleation
         387     4467528 :    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     4467528 :    real(r8), pointer :: t(:,:)          ! input temperature (K)
         396     4467528 :    real(r8), pointer :: qn(:,:)         ! input water vapor mixing ratio (kg/kg)
         397     4467528 :    real(r8), pointer :: qc(:,:)         ! cloud water mixing ratio (kg/kg)
         398     4467528 :    real(r8), pointer :: qi(:,:)         ! cloud ice mixing ratio (kg/kg)
         399     4467528 :    real(r8), pointer :: ni(:,:)         ! cloud ice number conc (1/kg)
         400     4467528 :    real(r8), pointer :: pmid(:,:)       ! pressure at layer midpoints (pa)
         401             :
         402     4467528 :    real(r8), pointer :: aer_mmr(:,:)    ! aerosol mass mixing ratio
         403             :
         404     4467528 :    real(r8), pointer :: ast(:,:)
         405             :    real(r8) :: icecldf(pcols,pver)  ! ice cloud fraction
         406     4467528 :    real(r8), pointer :: qsatfac(:,:)      ! Subgrid cloud water saturation scaling factor.
         407             :
         408             :    real(r8) :: rho(pcols,pver)      ! air density (kg m-3)
         409             :
         410     4467528 :    real(r8), allocatable :: naer2(:,:,:)    ! bulk aerosol number concentration (1/m3)
         411     4467528 :    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     4467528 :    real(r8), pointer :: amb_num(:,:)
         465     4467528 :    real(r8), pointer :: amb_mmr(:,:)
         466     4467528 :    real(r8), pointer :: cld_num(:,:)
         467     4467528 :    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     4467528 :    lchnk = state%lchnk
         477     4467528 :    ncol  = state%ncol
         478     4467528 :    t     => state%t
         479     4467528 :    qn    => state%q(:,:,1)
         480     4467528 :    qc    => state%q(:,:,cldliq_idx)
         481     4467528 :    qi    => state%q(:,:,cldice_idx)
         482     4467528 :    ni    => state%q(:,:,numice_idx)
         483     4467528 :    pmid  => state%pmid
         484             :
         485  6942019032 :    rho(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:))
         486             :
         487     4467528 :    if (clim_modal_aero) then
         488             :
         489     4467528 :       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     4467528 :    itim_old = pbuf_old_tim_idx()
         512    17870112 :    call pbuf_get_field(pbuf, ast_idx, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         513             :
         514  6942019032 :    icecldf(:ncol,:pver) = ast(:ncol,:pver)
         515             :
         516             :    ! naai and naai_hom are the outputs from this parameterization
         517     4467528 :    call pbuf_get_field(pbuf, naai_idx, naai)
         518     4467528 :    call pbuf_get_field(pbuf, naai_hom_idx, naai_hom)
         519  6942019032 :    naai(1:ncol,1:pver)     = 0._r8
         520  6942019032 :    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     4467528 :    troplev(:) = 0
         527             :    !REMOVECAM_END
         528     4467528 :    call tropopause_findChemTrop(state, troplev)
         529             :
         530     4467528 :    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     4467528 :    trop_pd(:,:) = 0._r8
         535             :
         536   379739880 :    do k = top_lev, pver
         537  6270643080 :       do i = 1, ncol
         538  5890903200 :          trop_pd(i, troplev(i)) = 1._r8
         539             :
         540  6266175552 :          if (k <= troplev(i)) then
         541  3033110807 :             if (nucleate_ice_subgrid_strat .eq. -1._r8) then
         542           0 :                subgrid(i, k) = 1._r8 / qsatfac(i, k)
         543             :             else
         544  3033110807 :                subgrid(i, k) = nucleate_ice_subgrid_strat
         545             :             end if
         546             :          else
         547  2857792393 :             if (nucleate_ice_subgrid .eq. -1._r8) then
         548           0 :                subgrid(i, k) = 1._r8 / qsatfac(i, k)
         549             :             else
         550  2857792393 :                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  6942019032 :    nihf(1:ncol,1:pver)  = 0._r8
         559  6942019032 :    niimm(1:ncol,1:pver) = 0._r8
         560  6942019032 :    nidep(1:ncol,1:pver) = 0._r8
         561  6942019032 :    nimey(1:ncol,1:pver) = 0._r8
         562             :
         563  6942019032 :    regm(1:ncol,1:pver) = 0._r8
         564             :
         565     4467528 :    if (use_preexisting_ice) then
         566     4467528 :       fhom(:,:)     = 0.0_r8
         567     4467528 :       wice(:,:)     = 0.0_r8
         568     4467528 :       weff(:,:)     = 0.0_r8
         569     4467528 :       INnso4(:,:)   = 0.0_r8
         570     4467528 :       INnbc(:,:)    = 0.0_r8
         571     4467528 :       INndust(:,:)  = 0.0_r8
         572     4467528 :       INondust(:,:)  = 0.0_r8
         573     4467528 :       INhet(:,:)    = 0.0_r8
         574     4467528 :       INhom(:,:)    = 0.0_r8
         575     4467528 :       INFrehom(:,:) = 0.0_r8
         576     4467528 :       INFreIN(:,:)  = 0.0_r8
         577             :    endif
         578             :
         579   379739880 :    do k = top_lev, pver
         580             :       ! Get humidity and saturation vapor pressures
         581   375272352 :       call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gammas(1:ncol))
         582             :
         583  6270643080 :       do i = 1, ncol
         584             :
         585  5890903200 :          relhum(i,k) = qn(i,k)/qs(i)
         586             :
         587             :          ! get cloud fraction, check for minimum
         588  6266175552 :          icldm(i,k) = max(icecldf(i,k), mincld)
         589             :
         590             :       end do
         591             :    end do
         592             :
         593     4467528 :    dust_num_col = 0._r8
         594     4467528 :    sulf_num_col = 0._r8
         595     4467528 :    sulf_num_tot_col = 0._r8
         596     4467528 :    soot_num_col = 0._r8
         597             :
         598     4467528 :    if (clim_modal_aero) then
         599             :
         600    13402584 :       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     4467528 :                                            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   379739880 :    kloop: do k = top_lev, pver
         619  6270643080 :       iloop: do i = 1, ncol
         620             :
         621  5890903200 :          so4_num_st_cr_tot = 0._r8
         622             :
         623  6266175552 :          freezing: if (t(i,k) < tmelt - 5._r8) then
         624             :
         625             :             ! set aerosol number for so4, soot, and dust with units #/cm^3
         626  4575144061 :             so4_num = sulf_num_col(i,k)
         627  4575144061 :             dst_num = dust_num_col(i,k)
         628  4575144061 :             so4_num_st_cr_tot=sulf_num_tot_col(i,k)
         629             :
         630             :             ! *** Turn off soot nucleation ***
         631  4575144061 :             soot_num = 0.0_r8
         632             :
         633  4575144061 :             if (cam_physpkg_is("cam7")) then
         634             :
         635             :                call nucleati( &
         636  9150288122 :                     wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k),   &
         637  4575144061 :                     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 13725432183 :                     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  4575144061 :             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 22875720305 :                do m = 1, aero_props%nbins()
         667             :
         668 22875720305 :                   if (aero_props%icenuc_updates_num(m)) then
         669             :
         670             :                      ! constituents of this bin will need to be updated
         671             :
         672  4575144061 :                      call aero_state%get_ambient_num(m, amb_num)
         673  4575144061 :                      call aero_state%get_cldbrne_num(m, cld_num)
         674             :
         675  4575144061 :                      if (amb_num(i,k)>0._r8) then
         676  4575144061 :                         delmmr_sum = 0._r8
         677  4575144061 :                         delnum_sum = 0._r8
         678             :
         679             :                         ! iterate over the species within the bin
         680 18300576244 :                         do l = 1, aero_props%nspecies(m)
         681 18300576244 :                            if (aero_props%icenuc_updates_mmr(m,l)) then
         682             :
         683  4575144061 :                               call aero_props%species_type(m, l, spectype)
         684  4575144061 :                               call aero_state%icenuc_size_wght( m, i,k, spectype, use_preexisting_ice, wght)
         685             :
         686  4575144061 :                               if (wght>0._r8) then
         687             :
         688             :                                  ! this aerosol constituent will be updated
         689             :
         690  4575144061 :                                  idxtmp = aer_cnst_idx(m,l)
         691             :
         692  4575144061 :                                  call aero_state%get_ambient_mmr(l,m,amb_mmr)
         693  4575144061 :                                  call aero_state%get_cldbrne_mmr(l,m,cld_mmr)
         694             :
         695             :                                  ! determine change in aerosol mass
         696  4575144061 :                                  delmmr = 0._r8
         697  4575144061 :                                  delnum = 0._r8
         698  4575144061 :                                  if (trim(spectype)=='dust') then
         699  4575144061 :                                     if (dst_num>0._r8) then
         700  4575144061 :                                        delmmr = (odst_num / dst_num) * icldm(i,k) * amb_mmr(i,k) * wght
         701  4575144061 :                                        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  4575144061 :                                  if (idxtmp>0) then
         711             :                                     ! constituent tendency (for transported species)
         712  4575144061 :                                     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  4575144061 :                                  cld_mmr(i,k) = cld_mmr(i,k) + delmmr
         718             :
         719  4575144061 :                                  delmmr_sum = delmmr_sum + delmmr
         720  4575144061 :                                  delnum_sum = delnum_sum + delnum
         721             :                               end if
         722             :                            end if
         723             :                         end do
         724             :
         725  4575144061 :                         idxtmp = aer_cnst_idx(m,0)
         726             :
         727             :                         ! update aerosol state bin and tendency for grid box i,k
         728  4575144061 :                         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  4575144061 :             if (nucleate_ice_use_troplev .and. clim_modal_aero) then
         756  4575144061 :                if ((k < troplev(i)) .and. (nucleate_ice_strat > 0._r8) .and. (oso4_num > 0._r8)) then
         757        2560 :                   dso4_num = max(0._r8, (nucleate_ice_strat*so4_num_st_cr_tot - oso4_num) * 1e6_r8 / rho(i,k))
         758        2560 :                   naai(i,k) = naai(i,k) + dso4_num
         759        2560 :                   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  4575144061 :             if (cam_physpkg_is("cam7")) then
         775             :                !Updates for pumas v1.21+
         776             :
         777  4575144061 :                naai_hom(i,k) = nihf(i,k)/dtime
         778  4575144061 :                naai(i,k)= naai(i,k)/dtime
         779             :
         780             :                ! output activated ice (convert from #/kg -> #/m3/s)
         781  4575144061 :                nihf(i,k)     = nihf(i,k) *rho(i,k)/dtime
         782  4575144061 :                niimm(i,k)    = niimm(i,k)*rho(i,k)/dtime
         783  4575144061 :                nidep(i,k)    = nidep(i,k)*rho(i,k)/dtime
         784  4575144061 :                nimey(i,k)    = nimey(i,k)*rho(i,k)/dtime
         785             :
         786  4575144061 :                if (use_preexisting_ice) then
         787  4575144061 :                   INnso4(i,k) =so4_num*1e6_r8/dtime  ! (convert from #/cm3 -> #/m3/s)
         788  4575144061 :                   INnbc(i,k)  =soot_num*1e6_r8/dtime
         789  4575144061 :                   INndust(i,k)=dst_num*1e6_r8/dtime
         790  4575144061 :                   INondust(i,k)=odst_num*1e6_r8/dtime
         791  4575144061 :                   INFreIN(i,k)=1.0_r8          ! 1,ice nucleation occur
         792  4575144061 :                   INhet(i,k) = (niimm(i,k) + nidep(i,k))   ! #/m3/s, nimey not in cirrus
         793  4575144061 :                   INhom(i,k) = nihf(i,k)                 ! #/m3/s
         794  4575144061 :                   if (INhom(i,k).gt.1e3_r8)   then ! > 1/L
         795       39701 :                      INFrehom(i,k)=1.0_r8       ! 1, hom freezing occur
         796             :                   endif
         797             :
         798             :                   ! exclude  no ice nucleaton
         799  4575144061 :                   if ((INFrehom(i,k) < 0.5_r8) .and. (INhet(i,k) < 1.0_r8))   then
         800  4565036680 :                      INnso4(i,k) =0.0_r8
         801  4565036680 :                      INnbc(i,k)  =0.0_r8
         802  4565036680 :                      INndust(i,k)=0.0_r8
         803  4565036680 :                      INondust(i,k)=0.0_r8
         804  4565036680 :                      INFreIN(i,k)=0.0_r8
         805  4565036680 :                      INhet(i,k) = 0.0_r8
         806  4565036680 :                      INhom(i,k) = 0.0_r8
         807  4565036680 :                      INFrehom(i,k)=0.0_r8
         808  4565036680 :                      wice(i,k) = 0.0_r8
         809  4565036680 :                      weff(i,k) = 0.0_r8
         810  4565036680 :                      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     4467528 :    if (.not. clim_modal_aero) then
         858           0 :       deallocate( &
         859             :            naer2, &
         860           0 :            maerosol)
         861             :    end if
         862             :
         863     4467528 :    if (cam_physpkg_is("cam7")) then
         864             :       ! Updates for PUMAS v1.21+
         865     4467528 :       call outfld('NIHFTEN',   nihf, pcols, lchnk)
         866     4467528 :       call outfld('NIIMMTEN', niimm, pcols, lchnk)
         867     4467528 :       call outfld('NIDEPTEN', nidep, pcols, lchnk)
         868     4467528 :       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     4467528 :    call outfld('NIREGM', regm, pcols, lchnk)
         876     4467528 :    call outfld('NISUBGRID', subgrid, pcols, lchnk)
         877     4467528 :    call outfld('NITROP_PD', trop_pd, pcols, lchnk)
         878             :
         879     4467528 :    if (use_preexisting_ice) then
         880     4467528 :       call outfld( 'fhom' , fhom, pcols, lchnk)
         881     4467528 :       call outfld( 'WICE' , wice, pcols, lchnk)
         882     4467528 :       call outfld( 'WEFF' , weff, pcols, lchnk)
         883     4467528 :       if (cam_physpkg_is("cam7")) then
         884             :          ! Updates for PUMAS v1.21+
         885     4467528 :          call outfld('INnso4TEN',INnso4 , pcols,lchnk)
         886     4467528 :          call outfld('INnbcTEN',INnbc  , pcols,lchnk)
         887     4467528 :          call outfld('INndustTEN',INndust, pcols,lchnk)
         888     4467528 :          call outfld('INondustTEN',INondust, pcols,lchnk)
         889     4467528 :          call outfld('INhetTEN',INhet  , pcols,lchnk)
         890     4467528 :          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     4467528 :       call outfld('INFrehom',INFrehom,pcols,lchnk)
         900     4467528 :       call outfld('INFreIN ',INFreIN, pcols,lchnk)
         901             :    end if
         902             :
         903     8935056 : end subroutine nucleate_ice_cam_calc
         904             :
         905             : !================================================================================================
         906             :
         907             : end module nucleate_ice_cam

Generated by: LCOV version 1.14