diff --git a/schemes/chemistry/aerosol_optics.F90 b/schemes/chemistry/aerosol_optics.F90 new file mode 100644 index 00000000..3b104dbe --- /dev/null +++ b/schemes/chemistry/aerosol_optics.F90 @@ -0,0 +1,552 @@ +! Aerosol optical properties for RRTMGP radiation. +! Computes per-band SW/LW aerosol optics from aerosol model instances. +! +! Author: Haipeng Lin, NSF-NCAR/CGD/AMP, March 2026 +module aerosol_optics + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: aerosol_optics_init + public :: aerosol_optics_run + + ! Water refractive indices (read in init, used in run) + complex(kind_phys), allocatable :: crefwsw(:) + complex(kind_phys), allocatable :: crefwlw(:) + +contains + +!> \section arg_table_aerosol_optics_init Argument Table +!! \htmlinclude aerosol_optics_init.html + subroutine aerosol_optics_init(N_DIAG, active_calls, & + water_refindex_file, & + nswbands, nlwbands, & + num_bulk_aer, & + errmsg, errflg) + use ccpp_io_reader, only: create_netcdf_reader_t + use ccpp_io_reader, only: abstract_netcdf_reader_t + use aerosol_mmr_host, only: rad_aer_diag_init + use radiative_aerosol_definitions, only: aerlist_t + use radiative_aerosol_definitions, only: bulk_aerosol_list + + use aerosol_instances_mod, only: aerosol_instances_is_active + + integer, intent(in) :: N_DIAG + logical, intent(in) :: active_calls(0:) + ! n.b.: the lower bound of the active_calls array has to be set here. + ! in radiative_aerosol_definitions, active_calls is 0-indexed because + ! 0 is the climate list index. + ! The Fortran 90 standard (and higher) states assumed-shape arrays + ! lower bound default to 1 unless specified (ISO/IEC 1539 5.1.2.4.2) + + !type(aerlist_t), intent(in) :: bulk_aerosol_list(:) + ! does not work: errors with + ! parse_source.CCPPError: No ddt_lib or ddt aerlist_t not in ddt_lib + + character(len=512), intent(in) :: water_refindex_file + integer, intent(in) :: nswbands + integer, intent(in) :: nlwbands + + integer, intent(out) :: num_bulk_aer + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: i + + ! Locals for reading water refractive index + class(abstract_netcdf_reader_t), pointer :: file_reader + real(kind_phys), allocatable :: refr(:), refi(:) + + errmsg = '' + errflg = 0 + + ! Set number of bulk aerosol constituents in climate list + num_bulk_aer = bulk_aerosol_list(0)%numaerosols + + ! Register aerosol diagnostic history fields for active lists + do i = 1, N_DIAG + if (active_calls(i)) then + call rad_aer_diag_init(bulk_aerosol_list(i)) + end if + end do + + !------------------------------------------------- + ! Read water refractive index data + !------------------------------------------------- + if (len_trim(water_refindex_file) == 0 .or. & + trim(water_refindex_file) == 'UNSET' .or. & + trim(water_refindex_file) == 'UNSET_PATH') then + ! if modal or carma aerosol are active this is a fatal error + if(aerosol_instances_is_active('modal') .or. & + aerosol_instances_is_active('carma')) then + errmsg = 'water_refindex_file must be specified when modal or carma aerosol models are active.' + errflg = 1 + return + end if + else + file_reader => create_netcdf_reader_t() + + call file_reader%open_file(trim(water_refindex_file), errmsg, errflg) + if (errflg /= 0) return + + ! Read shortwave refractive indices + call file_reader%get_var('refindex_real_water_sw', refr, errmsg, errflg) + if (errflg /= 0) return + call file_reader%get_var('refindex_im_water_sw', refi, errmsg, errflg) + if (errflg /= 0) return + + if (size(refr) /= nswbands) then + errflg = 1 + errmsg = 'aerosol_optics_init: nswbands mismatch in water refindex file' + return + end if + + allocate(crefwsw(nswbands)) + crefwsw(:) = cmplx(refr(:), refi(:), kind=kind_phys) + + deallocate(refr, refi) + + ! Read longwave refractive indices + call file_reader%get_var('refindex_real_water_lw', refr, errmsg, errflg) + if (errflg /= 0) return + call file_reader%get_var('refindex_im_water_lw', refi, errmsg, errflg) + if (errflg /= 0) return + + if (size(refr) /= nlwbands) then + errflg = 1 + errmsg = 'aerosol_optics_init: nlwbands mismatch in water refindex file' + return + end if + + allocate(crefwlw(nlwbands)) + crefwlw(:) = cmplx(refr(:), refi(:), kind=kind_phys) + + deallocate(refr, refi) + + call file_reader%close_file(errmsg, errflg) + if (errflg /= 0) return + + deallocate(file_reader) + end if + + end subroutine aerosol_optics_init + +!> \section arg_table_aerosol_optics_run Argument Table +!! \htmlinclude aerosol_optics_run.html + subroutine aerosol_optics_run( & + ncol, pver, & + const_props, & + nswbands, nlwbands, & + top_lev, & + rga, & + idx_sw_diag, & + num_bulk_aer, & + relh, pdeldry, constituents, & + aer_tau, aer_tau_w, aer_tau_w_g, aer_lw_abs, & + dustaod, sulfaod, bcaod, pomaod, soaaod, ssltaod, & + aodabsbc, & + burdendust, burdenso4, burdenbc, burdenpom, burdensoa, burdenseasalt, & + ssavis, aodvis, & + odv_col_aod, & + errmsg, errflg) + + ! framework dependency for const_props + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + + ! dependency to get constituent index + use ccpp_const_utils, only: ccpp_const_get_idx + + ! portable optics core: + use aerosol_optics_core, only: aerosol_optics_sw_bin, aerosol_optics_lw_bin + + ! host-model dependency for aerosol objects: + use phys_prop, only: nrh ! # of relative humidity bins for table lookup + use aerosol_instances_mod, only: aerosol_instances_get_props, & + aerosol_instances_get_state, & + aerosol_instances_get_num_models + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + use phys_prop, only: ot_length + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + type(ccpp_constituent_prop_ptr_t), & + intent(in) :: const_props(:) ! ccpp constituent properties pointer + integer, intent(in) :: nswbands ! number of SW bands [count] + integer, intent(in) :: nlwbands ! number of LW bands [count] + integer, intent(in) :: top_lev ! top level for aerosol model [index] + real(kind_phys), intent(in) :: rga + integer, intent(in) :: idx_sw_diag ! index of SW diagnostic band [index] + integer, intent(in) :: num_bulk_aer ! number of bulk aerosol constituents in climate list [count] + real(kind_phys), intent(in) :: relh(:, :) ! relative humidity [fraction] + real(kind_phys), intent(in) :: pdeldry(:, :) ! dry air pressure thickness [Pa] + real(kind_phys), intent(in) :: constituents(:, :, :) + + ! Output arguments + real(kind_phys), intent(out) :: aer_tau(:, :, :) ! SW extinction optical depth [1] + real(kind_phys), intent(out) :: aer_tau_w(:, :, :) ! SW ssa * tau [1] + real(kind_phys), intent(out) :: aer_tau_w_g(:, :, :) ! SW asy * ssa * tau [1] + real(kind_phys), intent(out) :: aer_lw_abs(:, :, :) ! LW absorption optical depth [1] + + ! Per-species diagnostic outputs: + real(kind_phys), intent(out) :: dustaod(:) ! Dust AOD at vis [1] + real(kind_phys), intent(out) :: sulfaod(:) ! Sulfate AOD at vis [1] + real(kind_phys), intent(out) :: bcaod(:) ! BC AOD at vis [1] + real(kind_phys), intent(out) :: pomaod(:) ! POM AOD at vis [1] + real(kind_phys), intent(out) :: soaaod(:) ! SOA AOD at vis [1] + real(kind_phys), intent(out) :: ssltaod(:) ! Seasalt AOD at vis [1] + real(kind_phys), intent(out) :: aodabsbc(:) ! BC absorption OD at vis [1] + real(kind_phys), intent(out) :: burdendust(:) ! Dust column burden [kg m-2] + real(kind_phys), intent(out) :: burdenso4(:) ! Sulfate column burden [kg m-2] + real(kind_phys), intent(out) :: burdenbc(:) ! BC column burden [kg m-2] + real(kind_phys), intent(out) :: burdenpom(:) ! POM column burden [kg m-2] + real(kind_phys), intent(out) :: burdensoa(:) ! SOA column burden [kg m-2] + real(kind_phys), intent(out) :: burdenseasalt(:) ! Seasalt column burden [kg m-2] + real(kind_phys), intent(out) :: ssavis(:) ! Column SSA at vis [1] + real(kind_phys), intent(out) :: aodvis(:) ! Column AOD at vis [1] + + ! Per-constituent column visible OD for bulk aerosol (ncol, num_bulk_aer) + real(kind_phys), intent(out) :: odv_col_aod(:,:) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: iaermod, ibin, nbins, iwav, ilev, icol + integer :: num_aero_models, numrh + integer :: bam_cnt + real(kind_phys) :: mass(ncol, pver) + + ! Index into constituents array for VOLC_RAD_GEOM + integer :: idx_volc_rad_geom + + ! Aerosol object pointers + class(aerosol_properties), pointer :: aeroprops + class(aerosol_state), pointer :: aerostate + + real(kind_phys) :: sulfwtpct(ncol, pver) + + ! Per-bin SW optics + real(kind_phys) :: tau_bin(ncol, pver, nswbands) + real(kind_phys) :: ssa_bin(ncol, pver, nswbands) + real(kind_phys) :: asm_bin(ncol, pver, nswbands) + real(kind_phys) :: pabs_vis(ncol, pver) + real(kind_phys) :: dopaer0_vis(ncol, pver) + + ! Per-bin LW optics + real(kind_phys) :: tau_lw_bin(ncol, pver, nlwbands) + real(kind_phys) :: absorp_bin(ncol, pver, nlwbands) + + ! Volcanic geometric radius (pointer into constituents) + real(kind_phys), target :: geometric_radius(ncol, pver) + logical :: has_volc_rad_geom + + character(len=ot_length) :: opticstype + + ! Per-species diagnostic locals + integer :: ispec + real(kind_phys) :: wetvol_d(ncol, pver), watervol_d(ncol, pver) + real(kind_phys) :: vol_d(ncol) + real(kind_phys) :: specdens_d, hygro_aer_d + character(len=32) :: spectype_d + real(kind_phys), pointer :: specmmr_d(:, :) + complex(kind_phys), pointer :: specrefindex_d(:) + real(kind_phys) :: scatdust_d(ncol), absdust_d(ncol), hygrodust_d(ncol) + real(kind_phys) :: scatbc_d(ncol), absbc_d(ncol), hygrobc_d(ncol) + real(kind_phys) :: scatsulf_d(ncol), abssulf_d(ncol), hygrosulf_d(ncol) + real(kind_phys) :: scatpom_d(ncol), abspom_d(ncol), hygropom_d(ncol) + real(kind_phys) :: scatsoa_d(ncol), abssoa_d(ncol), hygrosoa_d(ncol) + real(kind_phys) :: scatsslt_d(ncol), abssslt_d(ncol), hygrosslt_d(ncol) + real(kind_phys) :: scath2o_d, absh2o_d, sumscat_d, sumabs_d, sumhygro_d + real(kind_phys) :: dopaer_d, palb_d, aodc_d + + errmsg = '' + errflg = 0 + + ! Initialize outputs + aer_tau = 0._kind_phys + aer_tau_w = 0._kind_phys + aer_tau_w_g = 0._kind_phys + aer_lw_abs = 0._kind_phys + + ! Initialize per-constituent visible OD output for BAM + bam_cnt = 0 + odv_col_aod(:,:) = 0._kind_phys + + ! Initialize per-species diagnostic outputs + dustaod(:ncol) = 0._kind_phys + sulfaod(:ncol) = 0._kind_phys + bcaod(:ncol) = 0._kind_phys + pomaod(:ncol) = 0._kind_phys + soaaod(:ncol) = 0._kind_phys + ssltaod(:ncol) = 0._kind_phys + aodabsbc(:ncol) = 0._kind_phys + burdendust(:ncol) = 0._kind_phys + burdenso4(:ncol) = 0._kind_phys + burdenbc(:ncol) = 0._kind_phys + burdenpom(:ncol) = 0._kind_phys + burdensoa(:ncol) = 0._kind_phys + burdenseasalt(:ncol) = 0._kind_phys + ssavis(:ncol) = 0._kind_phys + aodvis(:ncol) = 0._kind_phys + + ! Layer mass [kg/m2] + mass(:ncol, :) = pdeldry(:ncol, :) * rga + + ! Check if VOLC_RAD_GEOM is available + idx_volc_rad_geom = -1 + call ccpp_const_get_idx(const_props, & + 'VOLC_RAD_GEOM', idx_volc_rad_geom, errmsg, errflg) + if(errflg /= 0) return + + has_volc_rad_geom = (idx_volc_rad_geom > 0) + if (has_volc_rad_geom) then + geometric_radius(:ncol, :pver) = constituents(:ncol, :pver, idx_volc_rad_geom) + end if + + num_aero_models = aerosol_instances_get_num_models() + + ! Loop over aerosol models in the climate list. + aerosol_model_loop: do iaermod = 1, num_aero_models + + aeroprops => aerosol_instances_get_props(iaermod, list_idx=0) + aerostate => aerosol_instances_get_state(iaermod, list_idx=0) + + if (.not. associated(aeroprops)) cycle + if (.not. associated(aerostate)) cycle + + nbins = aeroprops%nbins() + + ! Sulfate weight percent for hygroscopic_wtp optics + sulfwtpct(:ncol, :pver) = aerostate%wgtpct(ncol, pver) + + ! Loop over aerosol bins/species + bin_loop: do ibin = 1, nbins + + ! Check if volcanic_radius optics needs geometric_radius + call aeroprops%optics_params(bin_ndx=ibin, opticstype=opticstype) + if (index(opticstype, 'volcanic_radius') > 0 .and. .not. has_volc_rad_geom) then + errflg = 1 + errmsg = 'VOLC_RAD_GEOM constituent required for volcanic_radius optics but not found' + return + end if + + !------------------------------------------------- + ! Call the portable core for shortwave optics calculation per-bin: + !------------------------------------------------- + call aerosol_optics_sw_bin(aeroprops, aerostate, ibin, & + ncol, pver, top_lev, nswbands, nlwbands, nrh, & + idx_sw_diag, & + relh, sulfwtpct, mass, crefwsw, crefwlw, & + geometric_radius, & + tau_bin, ssa_bin, asm_bin, & + pabs_vis, dopaer0_vis, & + errmsg, errflg) + if (errflg /= 0) return + + ! Accumulate weighted products + wavelength: do iwav = 1, nswbands + vertical: do ilev = top_lev, pver + column: do icol = 1, ncol + aer_tau(icol, ilev, iwav) = aer_tau(icol, ilev, iwav) & + + tau_bin(icol, ilev, iwav) + aer_tau_w(icol, ilev, iwav) = aer_tau_w(icol, ilev, iwav) & + + tau_bin(icol, ilev, iwav) * ssa_bin(icol, ilev, iwav) + aer_tau_w_g(icol, ilev, iwav) = aer_tau_w_g(icol, ilev, iwav) & + + tau_bin(icol, ilev, iwav) * ssa_bin(icol, ilev, iwav) & + * asm_bin(icol, ilev, iwav) + end do column + end do vertical + end do wavelength + + !------------------------------------------------- + ! Per-species diagnostic accumulation at vis band + !------------------------------------------------- + wetvol_d(:ncol, :) = aerostate%wet_volume(aeroprops, ibin, ncol, pver) + watervol_d(:ncol, :) = aerostate%water_volume(aeroprops, ibin, ncol, pver) + + do ilev = top_lev, pver + ! Reset per-species accumulators for this bin/level + scatdust_d(:ncol) = 0._kind_phys; absdust_d(:ncol) = 0._kind_phys; hygrodust_d(:ncol) = 0._kind_phys + scatbc_d(:ncol) = 0._kind_phys; absbc_d(:ncol) = 0._kind_phys; hygrobc_d(:ncol) = 0._kind_phys + scatsulf_d(:ncol) = 0._kind_phys; abssulf_d(:ncol) = 0._kind_phys; hygrosulf_d(:ncol) = 0._kind_phys + scatpom_d(:ncol) = 0._kind_phys; abspom_d(:ncol) = 0._kind_phys; hygropom_d(:ncol) = 0._kind_phys + scatsoa_d(:ncol) = 0._kind_phys; abssoa_d(:ncol) = 0._kind_phys; hygrosoa_d(:ncol) = 0._kind_phys + scatsslt_d(:ncol) = 0._kind_phys; abssslt_d(:ncol) = 0._kind_phys; hygrosslt_d(:ncol) = 0._kind_phys + + do ispec = 1, aeroprops%nspecies(ibin) + call aeroprops%get(ibin, ispec, density=specdens_d, & + spectype=spectype_d, refindex_sw=specrefindex_d, hygro=hygro_aer_d) + call aerostate%get_ambient_mmr(species_ndx=ispec, bin_ndx=ibin, mmr=specmmr_d) + + do icol = 1, ncol + vol_d(icol) = specmmr_d(icol, ilev) / specdens_d + + select case (trim(spectype_d)) + case('dust') + burdendust(icol) = burdendust(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatdust_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + absdust_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygrodust_d(icol) = vol_d(icol)*hygro_aer_d + case('black-c') + burdenbc(icol) = burdenbc(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatbc_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + absbc_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygrobc_d(icol) = vol_d(icol)*hygro_aer_d + case('sulfate') + burdenso4(icol) = burdenso4(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatsulf_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + abssulf_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygrosulf_d(icol) = vol_d(icol)*hygro_aer_d + case('p-organic') + burdenpom(icol) = burdenpom(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatpom_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + abspom_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygropom_d(icol) = vol_d(icol)*hygro_aer_d + case('s-organic') + burdensoa(icol) = burdensoa(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatsoa_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + abssoa_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygrosoa_d(icol) = vol_d(icol)*hygro_aer_d + case('seasalt') + burdenseasalt(icol) = burdenseasalt(icol) + specmmr_d(icol, ilev)*mass(icol, ilev) + if (associated(specrefindex_d)) then + scatsslt_d(icol) = vol_d(icol)*specrefindex_d(idx_sw_diag)%re + abssslt_d(icol) = -vol_d(icol)*specrefindex_d(idx_sw_diag)%im + end if + hygrosslt_d(icol) = vol_d(icol)*hygro_aer_d + end select + end do ! icol + end do ! ispec + + ! Volume-based partitioning and per-species AOD accumulation + + ! FIXME: vol_d holds last species only (matches CAM behavior); + ! should be fixed to accumulate over all species + ! a corresponding change has to be made in aerosol_optics_cam.F90 in CAM. + ! This will all be revisited in the future in CAM as we move to runtime + ! configurability of which species are included in the optics calculation. + ! (hplin, 5/18/26) + do icol = 1, ncol + if (wetvol_d(icol, ilev) > 1.e-40_kind_phys .and. vol_d(icol) > 0._kind_phys) then + scath2o_d = watervol_d(icol, ilev)*crefwsw(idx_sw_diag)%re + absh2o_d = -watervol_d(icol, ilev)*crefwsw(idx_sw_diag)%im + sumscat_d = scatsulf_d(icol) + scatpom_d(icol) + scatsoa_d(icol) + scatbc_d(icol) + & + scatdust_d(icol) + scatsslt_d(icol) + scath2o_d + sumabs_d = abssulf_d(icol) + abspom_d(icol) + abssoa_d(icol) + absbc_d(icol) + & + absdust_d(icol) + abssslt_d(icol) + absh2o_d + sumhygro_d = hygrosulf_d(icol) + hygropom_d(icol) + hygrosoa_d(icol) + hygrobc_d(icol) + & + hygrodust_d(icol) + hygrosslt_d(icol) + + ! Partition water contribution proportionally to hygroscopicity + if (sumscat_d>0._kind_phys .and. sumabs_d>0._kind_phys .and. sumhygro_d>0._kind_phys) then + scatdust_d(icol) = (scatdust_d(icol) + scath2o_d*hygrodust_d(icol)/sumhygro_d)/sumscat_d + absdust_d(icol) = (absdust_d(icol) + absh2o_d*hygrodust_d(icol)/sumhygro_d)/sumabs_d + + scatsulf_d(icol) = (scatsulf_d(icol) + scath2o_d*hygrosulf_d(icol)/sumhygro_d)/sumscat_d + abssulf_d(icol) = (abssulf_d(icol) + absh2o_d*hygrosulf_d(icol)/sumhygro_d)/sumabs_d + + scatpom_d(icol) = (scatpom_d(icol) + scath2o_d*hygropom_d(icol)/sumhygro_d)/sumscat_d + abspom_d(icol) = (abspom_d(icol) + absh2o_d*hygropom_d(icol)/sumhygro_d)/sumabs_d + + scatsoa_d(icol) = (scatsoa_d(icol) + scath2o_d*hygrosoa_d(icol)/sumhygro_d)/sumscat_d + abssoa_d(icol) = (abssoa_d(icol) + absh2o_d*hygrosoa_d(icol)/sumhygro_d)/sumabs_d + + scatbc_d(icol) = (scatbc_d(icol) + scath2o_d*hygrobc_d(icol)/sumhygro_d)/sumscat_d + absbc_d(icol) = (absbc_d(icol) + absh2o_d*hygrobc_d(icol)/sumhygro_d)/sumabs_d + + scatsslt_d(icol) = (scatsslt_d(icol) + scath2o_d*hygrosslt_d(icol)/sumhygro_d)/sumscat_d + abssslt_d(icol) = (abssslt_d(icol) + absh2o_d*hygrosslt_d(icol)/sumhygro_d)/sumabs_d + end if + + dopaer_d = tau_bin(icol, ilev, idx_sw_diag) + palb_d = ssa_bin(icol, ilev, idx_sw_diag) + + ! Accumulate per-species AODs + aodabsbc(icol) = aodabsbc(icol) + absbc_d(icol)*dopaer_d*(1.0_kind_phys - palb_d) + + aodc_d = (abssulf_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatsulf_d(icol))*dopaer_d + sulfaod(icol) = sulfaod(icol) + aodc_d + + aodc_d = (abspom_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatpom_d(icol))*dopaer_d + pomaod(icol) = pomaod(icol) + aodc_d + + aodc_d = (abssoa_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatsoa_d(icol))*dopaer_d + soaaod(icol) = soaaod(icol) + aodc_d + + aodc_d = (absbc_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatbc_d(icol))*dopaer_d + bcaod(icol) = bcaod(icol) + aodc_d + + aodc_d = (abssslt_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatsslt_d(icol))*dopaer_d + ssltaod(icol) = ssltaod(icol) + aodc_d + + aodc_d = (absdust_d(icol)*(1.0_kind_phys - palb_d) + palb_d*scatdust_d(icol))*dopaer_d + dustaod(icol) = dustaod(icol) + aodc_d + end if + end do ! icol + end do ! ilev (per-species diagnostics) + + !------------------------------------------------- + ! Per-bin column visible tau for bulk aerosol (ODV_ diagnostics) + ! (was in aer_vis_diag_mod in CAM) + !------------------------------------------------- + if (aeroprops%model_is('BAM')) then + bam_cnt = bam_cnt + 1 + if (bam_cnt <= num_bulk_aer) then + ! TODO: Use hist_fld_active to gate this when available in CAM-SIMA + do icol = 1, ncol + odv_col_aod(icol, bam_cnt) = sum(tau_bin(icol, top_lev:pver, idx_sw_diag)) + end do + end if + end if + + !------------------------------------------------- + ! Call the portable core for longwave optics calculation per-bin: + !------------------------------------------------- + call aerosol_optics_lw_bin(aeroprops, aerostate, ibin, & + ncol, pver, nswbands, nlwbands, nrh, & + relh, sulfwtpct, mass, crefwsw, crefwlw, & + geometric_radius, & + tau_lw_bin, absorp_bin, & + errmsg, errflg) + if (errflg /= 0) return + + ! Accumulate LW absorption optical depth + do iwav = 1, nlwbands + do ilev = 1, pver + do icol = 1, ncol + aer_lw_abs(icol, ilev, iwav) = aer_lw_abs(icol, ilev, iwav) & + + tau_lw_bin(icol, ilev, iwav) + end do + end do + end do + + end do bin_loop ! ibin + end do aerosol_model_loop ! iaermod + + ! Compute column-level SSA and AOD at vis from total aer_tau/aer_tau_w + do icol = 1, ncol + aodvis(icol) = sum(aer_tau(icol, :, idx_sw_diag)) + if (aodvis(icol) > 1.e-10_kind_phys) then + ssavis(icol) = sum(aer_tau_w(icol, :, idx_sw_diag)) / aodvis(icol) + else + ! magic number from CAM aerosol_optics_cam.F90. + ssavis(icol) = 0.925_kind_phys + end if + end do + + end subroutine aerosol_optics_run + +end module aerosol_optics diff --git a/schemes/chemistry/aerosol_optics.meta b/schemes/chemistry/aerosol_optics.meta new file mode 100644 index 00000000..ceb5e587 --- /dev/null +++ b/schemes/chemistry/aerosol_optics.meta @@ -0,0 +1,263 @@ +[ccpp-table-properties] + name = aerosol_optics + type = scheme + +[ccpp-arg-table] + name = aerosol_optics_init + type = scheme +[ N_DIAG ] + standard_name = number_of_radiative_aerosol_diagnostic_lists + units = count + type = integer + dimensions = () + intent = in +[ active_calls ] + standard_name = flag_for_active_radiative_aerosol_diagnostic_list + units = flag + type = logical + dimensions = (index_of_climate_radiative_aerosol_list:number_of_radiative_aerosol_diagnostic_lists) + intent = in +[ water_refindex_file ] + standard_name = filename_of_water_refractive_index_data + units = none + type = character | kind = len=512 + dimensions = () + intent = in +[ nswbands ] + standard_name = number_of_bands_for_shortwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ nlwbands ] + standard_name = number_of_bands_for_longwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ num_bulk_aer ] + standard_name = number_of_bulk_aerosol_constituents_in_climate_radiative_aerosol_list + units = count + type = integer + dimensions = () + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = aerosol_optics_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ const_props ] + standard_name = ccpp_constituent_properties + units = none + type = ccpp_constituent_prop_ptr_t + dimensions = (number_of_ccpp_constituents) + intent = in +[ nswbands ] + standard_name = number_of_bands_for_shortwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ nlwbands ] + standard_name = number_of_bands_for_longwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = index_of_air_pressure_at_top_of_aerosol_model + units = index + type = integer + dimensions = () + intent = in +[ rga ] + standard_name = reciprocal_of_gravitational_acceleration + units = s2 m-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ idx_sw_diag ] + standard_name = index_of_shortwave_band + units = index + type = integer + dimensions = () + intent = in +[ num_bulk_aer ] + standard_name = number_of_bulk_aerosol_constituents_in_climate_radiative_aerosol_list + units = count + type = integer + dimensions = () + intent = in +[ relh ] + standard_name = relative_humidity + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pdeldry ] + standard_name = air_pressure_thickness_of_dry_air + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ constituents ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = in +[ aer_tau ] + standard_name = shortwave_aerosol_extinction_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_tau_w ] + standard_name = shortwave_aerosol_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_tau_w_g ] + standard_name = shortwave_aerosol_asymmetry_parameter_times_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_lw_abs ] + standard_name = longwave_aerosol_absorption_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_longwave_radiation) + intent = out +[ dustaod ] + standard_name = aerosol_optical_depth_from_dust + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ sulfaod ] + standard_name = aerosol_optical_depth_from_sulfate + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ bcaod ] + standard_name = aerosol_optical_depth_from_black_carbon + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ pomaod ] + standard_name = aerosol_optical_depth_from_primary_organic_matter + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ soaaod ] + standard_name = aerosol_optical_depth_from_secondary_organic_aerosol + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ssltaod ] + standard_name = aerosol_optical_depth_from_seasalt + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ aodabsbc ] + standard_name = absorption_aerosol_optical_depth_from_black_carbon + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdendust ] + standard_name = column_burden_from_dust + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdenso4 ] + standard_name = column_burden_from_sulfate + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdenbc ] + standard_name = column_burden_from_black_carbon + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdenpom ] + standard_name = column_burden_from_primary_organic_matter + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdensoa ] + standard_name = column_burden_from_secondary_organic_aerosol + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ burdenseasalt ] + standard_name = column_burden_from_seasalt + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ssavis ] + standard_name = aerosol_single_scatter_albedo_at_vis + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ aodvis ] + standard_name = aerosol_optical_depth_at_vis + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ odv_col_aod ] + standard_name = column_visible_optical_depth_per_bulk_aerosol + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, number_of_bulk_aerosol_constituents_in_climate_radiative_aerosol_list) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/chemistry/aerosol_optics_namelist.xml b/schemes/chemistry/aerosol_optics_namelist.xml new file mode 100644 index 00000000..a6af954a --- /dev/null +++ b/schemes/chemistry/aerosol_optics_namelist.xml @@ -0,0 +1,93 @@ + + + + + + + + + + char*512 + radiation + aerosol_optics_nl + filename_of_water_refractive_index_data + none + + Full pathname of dataset for water refractive indices used in aerosol optics. + + + ${DIN_LOC_ROOT}/atm/cam/physprops/water_refindex_rrtmg_c080910.nc + + + + diff --git a/schemes/chemistry/aerosol_optics_zero_stub.F90 b/schemes/chemistry/aerosol_optics_zero_stub.F90 new file mode 100644 index 00000000..bc4798f3 --- /dev/null +++ b/schemes/chemistry/aerosol_optics_zero_stub.F90 @@ -0,0 +1,40 @@ +! Zero stub for aerosol optical properties. +! Sets all aerosol optics arrays to zero (no aerosol), for use in +! test SDFs where full aerosol optics calculations are not needed. +module aerosol_optics_zero_stub + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: aerosol_optics_zero_stub_run + +contains + +!> \section arg_table_aerosol_optics_zero_stub_run Argument Table +!! \htmlinclude aerosol_optics_zero_stub_run.html +!! + subroutine aerosol_optics_zero_stub_run( & + aer_tau, aer_tau_w, aer_tau_w_g, aer_lw_abs, & + errmsg, errflg) + + ! Outputs + real(kind_phys), intent(out) :: aer_tau(:,:,:) !< SW extinction OD + real(kind_phys), intent(out) :: aer_tau_w(:,:,:) !< SW ssa*tau + real(kind_phys), intent(out) :: aer_tau_w_g(:,:,:) !< SW asy*ssa*tau + real(kind_phys), intent(out) :: aer_lw_abs(:,:,:) !< LW absorption OD + + character(len=*), intent(out) :: errmsg !< CCPP error message + integer, intent(out) :: errflg !< CCPP error flag + + errmsg = '' + errflg = 0 + + aer_tau(:,:,:) = 0._kind_phys + aer_tau_w(:,:,:) = 0._kind_phys + aer_tau_w_g(:,:,:) = 0._kind_phys + aer_lw_abs(:,:,:) = 0._kind_phys + + end subroutine aerosol_optics_zero_stub_run + +end module aerosol_optics_zero_stub diff --git a/schemes/chemistry/aerosol_optics_zero_stub.meta b/schemes/chemistry/aerosol_optics_zero_stub.meta new file mode 100644 index 00000000..f8c2360a --- /dev/null +++ b/schemes/chemistry/aerosol_optics_zero_stub.meta @@ -0,0 +1,43 @@ +[ccpp-table-properties] + name = aerosol_optics_zero_stub + type = scheme + +[ccpp-arg-table] + name = aerosol_optics_zero_stub_run + type = scheme +[ aer_tau ] + standard_name = shortwave_aerosol_extinction_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_tau_w ] + standard_name = shortwave_aerosol_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_tau_w_g ] + standard_name = shortwave_aerosol_asymmetry_parameter_times_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = out +[ aer_lw_abs ] + standard_name = longwave_aerosol_absorption_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_longwave_radiation) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/chemistry/prescribed_aerosol_deposition_flux_namelist.xml b/schemes/chemistry/prescribed_aerosol_deposition_flux_namelist.xml index 03de316b..3dde5046 100644 --- a/schemes/chemistry/prescribed_aerosol_deposition_flux_namelist.xml +++ b/schemes/chemistry/prescribed_aerosol_deposition_flux_namelist.xml @@ -161,10 +161,10 @@ The cycle year of the prescribed aerosol deposition flux data if aerodep_flx_type is 'CYCLICAL'. Format: YYYY - Default: 2000 + Default: 0 - 2000 + 0 diff --git a/schemes/chemistry/prescribed_aerosols.F90 b/schemes/chemistry/prescribed_aerosols.F90 index ccb72c8c..025408ce 100644 --- a/schemes/chemistry/prescribed_aerosols.F90 +++ b/schemes/chemistry/prescribed_aerosols.F90 @@ -138,6 +138,7 @@ subroutine prescribed_aerosols_register( & if (prescribed_aero_file == 'UNSET' .or. & len_trim(prescribed_aero_file) == 0) then has_prescribed_aerosols = .false. + allocate(aerosol_constituents(0)) if (amIRoot) then write(iulog,*) subname//': No prescribed aerosols specified' end if @@ -179,6 +180,7 @@ subroutine prescribed_aerosols_register( & if(aero_count == 0) then has_prescribed_aerosols = .false. + allocate(aerosol_constituents(0)) return end if diff --git a/schemes/chemistry/prescribed_ozone.F90 b/schemes/chemistry/prescribed_ozone.F90 index ff0dd0e9..c27e8bd2 100644 --- a/schemes/chemistry/prescribed_ozone.F90 +++ b/schemes/chemistry/prescribed_ozone.F90 @@ -65,6 +65,7 @@ subroutine prescribed_ozone_register( & write(iulog,*) 'prescribed_ozone_register: ozone is prescribed in: ' // trim(filename) end if else + allocate(ozone_constituents(0)) return end if diff --git a/schemes/chemistry/prescribed_volcanic_aerosol.F90 b/schemes/chemistry/prescribed_volcanic_aerosol.F90 index 20e0a7e3..d9970f1f 100644 --- a/schemes/chemistry/prescribed_volcanic_aerosol.F90 +++ b/schemes/chemistry/prescribed_volcanic_aerosol.F90 @@ -75,6 +75,7 @@ subroutine prescribed_volcanic_aerosol_register( & ! Check if prescribed volcanic aerosols are enabled if (prescribed_volcaero_file == 'UNSET' .or. & len_trim(prescribed_volcaero_file) == 0) then + allocate(volcaero_constituents(0)) if (amIRoot) then write(iulog,*) subname//': No prescribed volcanic aerosols specified' end if diff --git a/schemes/chemistry/rebin_seasalt.F90 b/schemes/chemistry/rebin_seasalt.F90 new file mode 100644 index 00000000..3798dfb5 --- /dev/null +++ b/schemes/chemistry/rebin_seasalt.F90 @@ -0,0 +1,188 @@ +! Rebins 4 sea salt bins +! sslt1, sslt2, sslt3, sslt4 +! into 2 bins for radiation: +! SSLTA (accumulation mode) and SSLTC (coarse mode). +! +! Original author: Francis Vitt, 9 May 2008 +! CCPPized: Haipeng Lin, Mar 2026 +module rebin_seasalt + use ccpp_kinds, only: kind_phys + + implicit none + private + + ! public CCPP-compliant subroutines + public :: rebin_seasalt_register + public :: rebin_seasalt_init + public :: rebin_seasalt_run + + ! Fraction of total sea salt mass in coarse mode [1] + real(kind_phys), parameter :: wgt_sscm = 6.0_kind_phys / 7.0_kind_phys + + ! Constituent indices (set during init) + integer :: idx_sslt1, idx_sslt2, idx_sslt3, idx_sslt4 + integer :: idx_sslta, idx_ssltc + + ! Flag for whether sea salt constituents are available + logical :: has_sslt = .false. + +contains + + ! Register SSLTA and SSLTC as non-advected CCPP constituents. +!> \section arg_table_rebin_seasalt_register Argument Table +!! \htmlinclude rebin_seasalt_register.html + subroutine rebin_seasalt_register( & + seasalt_constituents, & + errmsg, errflg) + + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + + ! Output arguments: + type(ccpp_constituent_properties_t), allocatable, intent(out) :: seasalt_constituents(:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + character(len=*), parameter :: subname = 'rebin_seasalt_register' + + errmsg = '' + errflg = 0 + + allocate(seasalt_constituents(2), stat=errflg, errmsg=errmsg) + if (errflg /= 0) then + errmsg = subname // ': ' // trim(errmsg) + return + end if + + ! SSLTA: sea salt accumulation mode + call seasalt_constituents(1)%instantiate( & + std_name = 'SSLTA', & + diag_name = 'SSLTA', & + long_name = 'sea salt accumulation mode mixing ratio', & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + min_value = 1.e-36_kind_phys, & + advected = .false., & + water_species = .false., & + mixing_ratio_type = 'dry', & + errcode = errflg, & + errmsg = errmsg) + if (errflg /= 0) return + + ! SSLTC: sea salt coarse mode + call seasalt_constituents(2)%instantiate( & + std_name = 'SSLTC', & + diag_name = 'SSLTC', & + long_name = 'sea salt coarse mode mixing ratio', & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + min_value = 1.e-36_kind_phys, & + advected = .false., & + water_species = .false., & + mixing_ratio_type = 'dry', & + errcode = errflg, & + errmsg = errmsg) + if (errflg /= 0) return + + end subroutine rebin_seasalt_register + + ! Initialize constituent indices for sea salt bins and rebinned modes. +!> \section arg_table_rebin_seasalt_init Argument Table +!! \htmlinclude rebin_seasalt_init.html + subroutine rebin_seasalt_init( & + amIRoot, iulog, & + pcnst, const_props, & + errmsg, errflg) + + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_const_utils, only: ccpp_const_get_idx + + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + integer, intent(in) :: pcnst + type(ccpp_constituent_prop_ptr_t), & + intent(in) :: const_props(:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + character(len=*), parameter :: subname = 'rebin_seasalt_init' + + errmsg = '' + errflg = 0 + + ! Look up the 4 input sea salt bin constituent indices. + ! These are registered by prescribed_aerosols (bulk aerosol specifier). + call ccpp_const_get_idx(const_props, 'sslt1', idx_sslt1, errmsg, errflg) + if (errflg /= 0) return + call ccpp_const_get_idx(const_props, 'sslt2', idx_sslt2, errmsg, errflg) + if (errflg /= 0) return + call ccpp_const_get_idx(const_props, 'sslt3', idx_sslt3, errmsg, errflg) + if (errflg /= 0) return + call ccpp_const_get_idx(const_props, 'sslt4', idx_sslt4, errmsg, errflg) + if (errflg /= 0) return + + if(idx_sslt1 < 0 .or. idx_sslt2 < 0 .or. & + idx_sslt3 < 0 .or. idx_sslt4 < 0) then + if(amIRoot) then + write(iulog,*) subname//': did not find sslt1/2/3/4 and will not rebin seasalt to SSLTA/SSLTC' + end if + return + end if + + ! Look up the output rebinned constituent indices (registered by this scheme) + call ccpp_const_get_idx(const_props, & + 'SSLTA', idx_sslta, errmsg, errflg) + if (errflg /= 0) return + call ccpp_const_get_idx(const_props, & + 'SSLTC', idx_ssltc, errmsg, errflg) + if (errflg /= 0) return + + has_sslt = .true. + + if (amIRoot) then + write(iulog,*) subname//': initialized; coarse mode weight = ', wgt_sscm + end if + + end subroutine rebin_seasalt_init + + ! Rebin 4 sea salt bins into accumulation and coarse modes. +!> \section arg_table_rebin_seasalt_run Argument Table +!! \htmlinclude rebin_seasalt_run.html + subroutine rebin_seasalt_run( & + ncol, pver, pcnst, & + constituents, & + errmsg, errflg) + + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pcnst + real(kind_phys), intent(inout) :: constituents(:,:,:) ! constituent mixing ratios [kg kg-1] + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + real(kind_phys) :: sslt_sum ! sum of 4 sea salt bins [kg kg-1] + + errmsg = '' + errflg = 0 + + if (.not. has_sslt) return + + do k = 1, pver + do i = 1, ncol + sslt_sum = constituents(i, k, idx_sslt1) & + + constituents(i, k, idx_sslt2) & + + constituents(i, k, idx_sslt3) & + + constituents(i, k, idx_sslt4) + + ! Accumulation mode: fraction (1 - wgt_sscm) = 1/7 of total + constituents(i, k, idx_sslta) = (1.0_kind_phys - wgt_sscm) * sslt_sum + + ! Coarse mode: fraction wgt_sscm = 6/7 of total + constituents(i, k, idx_ssltc) = wgt_sscm * sslt_sum + end do + end do + + end subroutine rebin_seasalt_run + +end module rebin_seasalt diff --git a/schemes/chemistry/rebin_seasalt.meta b/schemes/chemistry/rebin_seasalt.meta new file mode 100644 index 00000000..348d8b2a --- /dev/null +++ b/schemes/chemistry/rebin_seasalt.meta @@ -0,0 +1,106 @@ +[ccpp-table-properties] + name = rebin_seasalt + type = scheme + +[ccpp-arg-table] + name = rebin_seasalt_register + type = scheme +[ seasalt_constituents ] + standard_name = dynamic_constituents_for_rebin_seasalt + units = none + type = ccpp_constituent_properties_t + dimensions = (:) + allocatable = True + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = rebin_seasalt_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ pcnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ const_props ] + standard_name = ccpp_constituent_properties + units = none + type = ccpp_constituent_prop_ptr_t + dimensions = (number_of_ccpp_constituents) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = rebin_seasalt_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pcnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ constituents ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = inout +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.F90 b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.F90 new file mode 100644 index 00000000..31ee96b1 --- /dev/null +++ b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.F90 @@ -0,0 +1,56 @@ +! Stub for heterogeneous freezing by classical nucleation theory (hetfrz_classnuc). +! +! This scheme exists to read the use_hetfrz_classnuc namelist flag and expose it +! via its CCPP standard name so downstream schemes (e.g. nucleate_ice_ccpp, +! eventually PUMAS) can see the flag. It performs no heterogeneous-freezing +! computation -- the real hetfrz_classnuc_ccpp is not yet ported to CAM-SIMA. +! +! Safe to use with use_hetfrz_classnuc=.true. for snapshot-testing flows where +! the flag must propagate downstream; a warning is printed at init time. +module hetfrz_classnuc_stub + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: hetfrz_classnuc_stub_init + +contains + +!> \section arg_table_hetfrz_classnuc_stub_init Argument Table +!! \htmlinclude hetfrz_classnuc_stub_init.html + subroutine hetfrz_classnuc_stub_init( & + amIRoot, iulog, & + use_hetfrz_classnuc, & + errmsg, errflg) + + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + logical, intent(in) :: use_hetfrz_classnuc + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + if (amIRoot) then + if (use_hetfrz_classnuc) then + write(iulog, '(A)') 'WARNING: hetfrz_classnuc_stub: use_hetfrz_classnuc=.true., but ' // & + 'heterogeneous freezing by classical nucleation theory is not yet ' // & + 'implemented in CAM-SIMA (stub scheme).' + write(iulog, '(A)') 'WARNING: hetfrz_classnuc_stub: The flag is still propagated to ' // & + 'downstream schemes (e.g. nucleate_ice_ccpp) via its CCPP standard ' // & + 'name, but no hetfrz_classnuc computation will occur.' + + ! this is not intended to error out for now as the flag needs to be passed to downstream schemes + ! (e.g., PUMAS) and modifies their behavior + else + write(iulog, '(A)') 'hetfrz_classnuc_stub: use_hetfrz_classnuc=.false. (stub scheme, ' // & + 'no-op).' + end if + end if + + end subroutine hetfrz_classnuc_stub_init + +end module hetfrz_classnuc_stub diff --git a/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.meta b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.meta new file mode 100644 index 00000000..fabd58a2 --- /dev/null +++ b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub.meta @@ -0,0 +1,37 @@ +[ccpp-table-properties] + name = hetfrz_classnuc_stub + type = scheme + +[ccpp-arg-table] + name = hetfrz_classnuc_stub_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ use_hetfrz_classnuc ] + standard_name = do_heterogeneous_freezing_by_classical_nucleation_theory + units = flag + type = logical + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/ccpp/hetfrz_classnuc_stub_namelist.xml b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub_namelist.xml new file mode 100644 index 00000000..35665f31 --- /dev/null +++ b/schemes/microp_aero/ccpp/hetfrz_classnuc_stub_namelist.xml @@ -0,0 +1,28 @@ + + + + + + + + logical + microphys + hetfrz_classnuc_nl + do_heterogeneous_freezing_by_classical_nucleation_theory + flag + + Switch to turn on heterogeneous freezing by classical nucleation theory. + Currently a stub in CAM-SIMA: the flag is propagated to downstream schemes + (nucleate_ice_ccpp, eventually PUMAS) but no hetfrz computation is performed. + Setting this to TRUE is supported for snapshot-testing flows where the flag + must match CAM; a warning is emitted at init time. + Default: FALSE; TRUE for CAM6 and CAM7. + + + .false. + .true. + .true. + + + + diff --git a/schemes/microp_aero/ccpp/ndrop_bam_ccpp.F90 b/schemes/microp_aero/ccpp/ndrop_bam_ccpp.F90 new file mode 100644 index 00000000..f0f9a60c --- /dev/null +++ b/schemes/microp_aero/ccpp/ndrop_bam_ccpp.F90 @@ -0,0 +1,163 @@ +! CCPP wrapper for BAM droplet activation (ndrop_bam). +! Looks up BAM aerosol objects from aerosol_instances and delegates to the +! portable ndrop_bam driver. +module ndrop_bam_ccpp + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: ndrop_bam_ccpp_init + public :: ndrop_bam_ccpp_run + + ! smallest mixing ratio considered in microphysics + real(kind_phys), parameter :: qsmall = 1.e-18_kind_phys + +contains + +!> \section arg_table_ndrop_bam_ccpp_init Argument Table +!! \htmlinclude ndrop_bam_ccpp_init.html + subroutine ndrop_bam_ccpp_init(& + amIRoot, iulog, & + mwh2o, r_universal, tmelt, rhoh2o, & + naer_all, psat, & + errmsg, errflg) + + use ndrop_bam, only: ndrop_bam_init + use ndrop_bam, only: psat_driver => psat + + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + real(kind_phys), intent(in) :: mwh2o + real(kind_phys), intent(in) :: r_universal + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: rhoh2o + integer, intent(out) :: naer_all + integer, intent(out) :: psat + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + call ndrop_bam_init(amIRoot, iulog, mwh2o, r_universal, tmelt, rhoh2o, & + naer_all) + + psat = psat_driver + + end subroutine ndrop_bam_ccpp_init + +!> \section arg_table_ndrop_bam_ccpp_run Argument Table +!! \htmlinclude ndrop_bam_ccpp_run.html + subroutine ndrop_bam_ccpp_run( & + ncol, pver, top_lev, & + gravit, rair, tmelt, cpair, rh2o, rhoh2o, latvap, & + pmid, t, wsub, qcld, ast, numliq, deltatin, & + npccn, nacon, ccn, naer2_diag, & + errmsg, errflg) + + use ndrop_bam, only: ndrop_bam_calc + use aerosol_instances_mod, only: aerosol_instances_get_props, & + aerosol_instances_get_state, & + aerosol_instances_get_num_models + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev + real(kind_phys), intent(in) :: gravit + real(kind_phys), intent(in) :: rair + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: cpair + real(kind_phys), intent(in) :: rh2o + real(kind_phys), intent(in) :: rhoh2o + real(kind_phys), intent(in) :: latvap + real(kind_phys), intent(in) :: pmid(:, :) + real(kind_phys), intent(in) :: t(:, :) + real(kind_phys), intent(in) :: wsub(:, :) + real(kind_phys), intent(in) :: qcld(:, :) + real(kind_phys), intent(in) :: ast(:, :) + real(kind_phys), intent(in) :: numliq(:, :) + real(kind_phys), intent(in) :: deltatin + real(kind_phys), intent(out) :: npccn(:, :) + real(kind_phys), intent(out) :: nacon(:, :, :) ! ncol, pver, num_dust=4 + real(kind_phys), intent(out) :: ccn(:, :, :) ! ncol, pver, psat + real(kind_phys), intent(out) :: naer2_diag(:, :, :) ! ncol, pver, naer_all + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k, iaermod + real(kind_phys) :: rho(ncol, pver) + class(aerosol_properties), pointer :: aero_props_bam + class(aerosol_state), pointer :: aero_state_bam + + errmsg = '' + errflg = 0 + + ! Find BAM properties and state from aerosol instances + aero_props_bam => null() + aero_state_bam => null() + do iaermod = 1, aerosol_instances_get_num_models() + aero_props_bam => aerosol_instances_get_props(iaermod, 0) + if (associated(aero_props_bam)) then + if (aero_props_bam%model_is('BAM')) then + aero_state_bam => aerosol_instances_get_state(iaermod, list_idx=0) + exit + end if + end if + aero_props_bam => null() + end do + + ! No BAM model found but microp_aero is active (i.e., aquaplanet). + ! Zero outputs and return. + if (.not. associated(aero_props_bam) .or. & + .not. associated(aero_state_bam)) then + npccn(:, :) = 0._kind_phys + nacon(:, :, :) = 0._kind_phys + ccn(:, :, :) = 0._kind_phys + naer2_diag(:, :, :) = 0._kind_phys + return + end if + + do k = top_lev, pver + do i = 1, ncol + rho(i,k) = pmid(i,k)/(rair*t(i,k)) + end do + end do + + call ndrop_bam_calc( & + aero_state = aero_state_bam, & + aero_props = aero_props_bam, & + ncol = ncol, & + pver = pver, & + top_lev = top_lev, & + gravit = gravit, & + rair = rair, & + tmelt = tmelt, & + cpair = cpair, & + rh2o = rh2o, & + rhoh2o = rhoh2o, & + latvap = latvap, & + rho = rho, & + tair = t, & + wsub = wsub, & + qcld = qcld, & + qsmall_in = qsmall, & + ast = ast, & + numliq = numliq, & + deltatin = deltatin, & ! below output: + npccn = npccn, & + nacon = nacon, & + ccn = ccn, & + naer2_diag = naer2_diag, & + errmsg = errmsg, & + errflg = errflg) + if(errflg /= 0) return + + end subroutine ndrop_bam_ccpp_run + +end module ndrop_bam_ccpp diff --git a/schemes/microp_aero/ccpp/ndrop_bam_ccpp.meta b/schemes/microp_aero/ccpp/ndrop_bam_ccpp.meta new file mode 100644 index 00000000..376a312d --- /dev/null +++ b/schemes/microp_aero/ccpp/ndrop_bam_ccpp.meta @@ -0,0 +1,216 @@ +[ccpp-table-properties] + name = ndrop_bam_ccpp + type = scheme + dependencies = ../ndrop_bam.F90 + +[ccpp-arg-table] + name = ndrop_bam_ccpp_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ mwh2o ] + standard_name = molecular_weight_of_h2o + units = g mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ r_universal ] + standard_name = universal_gas_constant + units = J K-1 kmol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhoh2o ] + standard_name = fresh_liquid_water_density_at_0c + units = kg m-3 + type = real | kind = kind_phys + dimensions = () + intent = in +[ naer_all ] + standard_name = number_of_bulk_aerosol_types + units = count + type = integer + dimensions = () + intent = out +[ psat ] + standard_name = number_of_ccn_supersaturation_levels + units = count + type = integer + dimensions = () + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = ndrop_bam_ccpp_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ cpair ] + standard_name = specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rh2o ] + standard_name = gas_constant_of_water_vapor + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhoh2o ] + standard_name = fresh_liquid_water_density_at_0c + units = kg m-3 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latvap ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qcld ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ ast ] + standard_name = stratiform_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ numliq ] + standard_name = mass_number_concentration_of_cloud_liquid_wrt_moist_air_and_condensed_water + units = kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ deltatin ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ npccn ] + # tendency of cloud droplet number from aerosol activation (#/kg/s) + standard_name = tendency_of_activated_cloud_condensation_nuclei_mass_number_concentration + units = kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ nacon ] + # contact freezing dust number concentration by bin (#/m3), 4 dust bins + standard_name = dust_number_concentration_by_size_bin_for_contact_freezing + units = m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, dust_size_bin_dimension) + intent = out +[ ccn ] + # CCN number concentration at fixed supersaturation levels (#/cm3), psat=6 levels + standard_name = cloud_condensation_nuclei_number_concentration_at_supersaturation + units = cm-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccn_supersaturation_levels) + intent = out +[ naer2_diag ] + # bulk aerosol number concentration diagnostic (#/m3), naer_all types + standard_name = bulk_aerosol_number_concentration_diagnostic + units = m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bulk_aerosol_types) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90 b/schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90 new file mode 100644 index 00000000..ed7f9fc6 --- /dev/null +++ b/schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90 @@ -0,0 +1,671 @@ +! Unified CCPP scheme for ice nucleation via the nucleati kernel. +! Uses the abstract interface to support BAM, MAM, and CARMA. +! +! This scheme replaces nucleate_ice_cam_calc for CAM-SIMA. +module nucleate_ice_ccpp + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: nucleate_ice_ccpp_init + public :: nucleate_ice_ccpp_run + + ! Module state from init + logical :: use_preexisting_ice_ = .false. + logical :: nucleate_ice_use_troplev_ = .false. + real(kind_phys) :: nucleate_ice_subgrid_ + real(kind_phys) :: nucleate_ice_subgrid_strat_ + real(kind_phys) :: nucleate_ice_strat_ + + ! Hardcoded scheme constants (previously CAM-side parameters/flags). + ! mincld_: floor for cloud fraction used in ice nucleation. CAM value was + ! real(r8), parameter :: mincld = 0.0001_r8 in microp_aero.F90. + ! use_nucleati_tendencies_: PUMAS v1.21+ (CAM7+) tendency path. Only supported + ! path in the CCPP version; legacy path will abort as not implemented. + real(kind_phys), parameter :: mincld_ = 0.0001_kind_phys + logical, parameter :: use_nucleati_tendencies_ = .true. + + ! Aerosol model flags (determined at init from aerosol_instances) + logical :: clim_modal_aero_ = .false. + logical :: clim_carma_aero_ = .false. + + ! index into aerosol_instances + ! NOTE: this replicates the microp_aero flow where only one aerosol + ! model/state is selected for use by microp_aero (ice nucleation, activation) + ! even if multiple aerosol models are active at the same time. + ! If modal or CARMA are active, they are used, otherwise bulk. (hplin, 4/7/26) + ! + ! This module also can run without an aerosol model, in aquaplanet compsets + ! with two-moment microphysics; nucleatei runs Meyers nucleation deposition + ! even without aerosol, in which case iaermod_selected_ is -1. + integer :: iaermod_selected_ = -1 + + ! Modal/CARMA-specific: whether cloud-borne tendencies are needed + logical :: clim_modal_carma_ = .false. + logical :: prog_modal_aero_ = .false. + + ! Aerosol bins/species dimensions (set in init, used in run for transfer) + integer :: nbins_ = 0 + integer :: nmaxspc_ = 0 + + ! Constituent tendency indices for interstitial-to-cloud-borne aerosol transfer. + ! Allocated in init if clim_modal_carma_ .and. use_preexisting_ice_. + ! aer_cnst_idx_(ibin, 0) = index for bin number/MMR (species 0) + ! aer_cnst_idx_(ibin, ispc) = index for species ispc MMR + ! Index is into CCPP constituent array; -1 means not transported. + integer, allocatable :: aer_cnst_idx_(:,:) + +contains + +!> \section arg_table_nucleate_ice_ccpp_init Argument Table +!! \htmlinclude nucleate_ice_ccpp_init.html + subroutine nucleate_ice_ccpp_init( & + iulog, & + const_props, & + use_preexisting_ice, use_hetfrz_classnuc, & + nucleate_ice_incloud, & + pi, & + nucleate_ice_subgrid, nucleate_ice_subgrid_strat, & + nucleate_ice_strat, nucleate_ice_use_troplev, & + prog_modal_aero, & + errmsg, errflg) + + use nucleate_ice, only: nucleati_init + use aerosol_instances_mod, only: aerosol_instances_is_active, & + aerosol_instances_get_props, & + aerosol_instances_get_num_models + use aerosol_properties_mod, only: aerosol_properties + + ! framework dependency for const_props + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_const_utils, only: ccpp_const_get_idx + + ! Input arguments + integer, intent(in) :: iulog + type(ccpp_constituent_prop_ptr_t), & + intent(in) :: const_props(:) + logical, intent(in) :: use_preexisting_ice + logical, intent(in) :: use_hetfrz_classnuc + logical, intent(in) :: nucleate_ice_incloud + real(kind_phys), intent(in) :: pi + real(kind_phys), intent(in) :: nucleate_ice_subgrid + real(kind_phys), intent(in) :: nucleate_ice_subgrid_strat + real(kind_phys), intent(in) :: nucleate_ice_strat + logical, intent(in) :: nucleate_ice_use_troplev + logical, intent(in) :: prog_modal_aero + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: iaermod + integer :: ibin, ispc, idxtmp + class(aerosol_properties), pointer :: aprops + character(len=32) :: tmpname + + errmsg = '' + errflg = 0 + + ! Store namelist parameters as module state + use_preexisting_ice_ = use_preexisting_ice + nucleate_ice_use_troplev_ = nucleate_ice_use_troplev + nucleate_ice_subgrid_ = nucleate_ice_subgrid + nucleate_ice_subgrid_strat_ = nucleate_ice_subgrid_strat + nucleate_ice_strat_ = nucleate_ice_strat + prog_modal_aero_ = prog_modal_aero + + ! Query aerosol model availability + clim_modal_aero_ = aerosol_instances_is_active('modal') + clim_carma_aero_ = aerosol_instances_is_active('carma') + + ! Select aerosol model (modal > CARMA > BAM > none) + iaermod_selected_ = -1 + do iaermod = 1, aerosol_instances_get_num_models() + aprops => aerosol_instances_get_props(iaermod, list_idx=0) + if (.not. associated(aprops)) cycle + + if (aprops%model_is('modal') .or. aprops%model_is('CARMA')) then + iaermod_selected_ = iaermod + clim_modal_carma_ = .true. + exit + else if (aprops%model_is('BAM')) then + iaermod_selected_ = iaermod + end if + end do + + ! Set up aerosol dimensions and constituent indices for transfer + nbins_ = 0 + nmaxspc_ = 0 + if (iaermod_selected_ > 0) then + aprops => aerosol_instances_get_props(iaermod_selected_, list_idx=0) + if (associated(aprops)) then + nbins_ = aprops%nbins() + nmaxspc_ = maxval(aprops%nspecies()) + end if + end if + + ! Set up constituent tendency indices for interstitial-to-cloud-borne + ! aerosol transfer (MAM/CARMA + use_preexisting_ice only). + ! Replaces nucleate_ice_cam_init L178-222. + if (clim_modal_carma_ .and. use_preexisting_ice_ .and. associated(aprops)) then + + allocate(aer_cnst_idx_(nbins_, 0:nmaxspc_), stat=errflg, errmsg=errmsg) + if(errflg /= 0) return + aer_cnst_idx_ = -1 + + do ibin = 1, nbins_ + if (aprops%icenuc_updates_num(ibin)) then + + ! Get constituent name for this bin (number or bin MMR) + if (aprops%icenuc_updates_mmr(ibin, 0)) then + call aprops%amb_mmr_name(ibin, 0, tmpname) + else + call aprops%amb_num_name(ibin, tmpname) + end if + + ! Look up CCPP constituent index by name (standard_name = constituent name) + call ccpp_const_get_idx(const_props, trim(tmpname), idxtmp, errmsg, errflg) + if (errflg /= 0) return + aer_cnst_idx_(ibin, 0) = idxtmp + + ! Iterate over species within the bin + do ispc = 1, aprops%nspecies(ibin) + if (aprops%icenuc_updates_mmr(ibin, ispc)) then + call aprops%amb_mmr_name(ibin, ispc, tmpname) + call ccpp_const_get_idx(const_props, trim(tmpname), idxtmp, errmsg, errflg) + if (errflg /= 0) return + aer_cnst_idx_(ibin, ispc) = idxtmp + end if + end do + + end if + end do + + end if + + ! Initialize the nucleati kernel + call nucleati_init(use_preexisting_ice, use_hetfrz_classnuc, & + nucleate_ice_incloud, iulog, pi, mincld_) + + end subroutine nucleate_ice_ccpp_init + +!> \section arg_table_nucleate_ice_ccpp_run Argument Table +!! \htmlinclude nucleate_ice_ccpp_run.html + subroutine nucleate_ice_ccpp_run( & + ncol, pver, pcnst, top_lev, dtime, & + t, pmid, qv, qc, qi, ni, & + rair, tmelt, pi, & + wsubi, & + aist, & + tropLev_chem, & + qsatfac, & + lat, & + naai, naai_hom, & + nihf, niimm, nidep, nimey, & + regm, subgrid_diag, trop_pd, & + fhom, wice, weff, & + INnso4, INnbc, INndust, INondust, & + INhet, INhom, INFrehom, INFreIN, & + ptend_q, & + errmsg, errflg) + + ! portable core science code: + use nucleate_ice, only: nucleati + + ! to_be_ccppized dependency: + use wv_saturation, only: qsat_water + + ! abstract aerosol interface: + use aerosol_instances_mod, only: aerosol_instances_get_props, & + aerosol_instances_get_state, & + aerosol_instances_get_num_models + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pcnst + integer, intent(in) :: top_lev ! top vertical level for cloud physics [index] + real(kind_phys), intent(in) :: dtime ! timestep [s] + real(kind_phys), intent(in) :: t(:, :) ! temperature [K] + real(kind_phys), intent(in) :: pmid(:, :) ! pressure at layer midpoints [Pa] + real(kind_phys), intent(in) :: qv(:, :) ! water vapor mixing ratio [kg kg-1] + real(kind_phys), intent(in) :: qc(:, :) ! cloud liquid mixing ratio [kg kg-1] + real(kind_phys), intent(in) :: qi(:, :) ! cloud ice mixing ratio [kg kg-1] + real(kind_phys), intent(in) :: ni(:, :) ! cloud ice number concentration [1 kg-1] + + real(kind_phys), intent(in) :: rair ! gas constant for dry air [J kg-1 K-1] + real(kind_phys), intent(in) :: tmelt ! freezing point of water [K] + real(kind_phys), intent(in) :: pi ! pi + + real(kind_phys), intent(in) :: wsubi(:, :) ! subgrid vertical velocity for ice nucleation [m s-1] + real(kind_phys), intent(in) :: aist(:, :) ! ice cloud fraction [fraction] + integer, intent(in) :: tropLev_chem(:)! chemical tropopause level [index] + real(kind_phys), intent(in) :: qsatfac(:, :) ! subgrid saturation scaling factor [1] + real(kind_phys), intent(in) :: lat(:) ! latitude [radians] + + ! Note: below units for naai...nimey are tendencies (with s-1) + ! as is convention for CAM7 PUMASv1.21+ + ! The legacy path is not supported in CCPP version of nucleate_ice; if it were + ! to be supported the units will be concentrations (no s-1); they would not be + ! divided by /dtime in this run phase. + + ! Output arguments + real(kind_phys), intent(out) :: naai(:, :) ! nucleated ice number [kg-1 s-1] + real(kind_phys), intent(out) :: naai_hom(:, :) ! nucleated ice number, hom. freezing only [kg-1 s-1] + + ! Diagnostic outputs + real(kind_phys), intent(out) :: nihf(:, :) ! ice nuclei from hom. freezing [m-3 s-1] + real(kind_phys), intent(out) :: niimm(:, :) ! ice nuclei from immersion freezing [m-3 s-1] + real(kind_phys), intent(out) :: nidep(:, :) ! ice nuclei from deposition nucleation [m-3 s-1] + real(kind_phys), intent(out) :: nimey(:, :) ! ice nuclei from Meyers deposition [m-3 s-1] + real(kind_phys), intent(out) :: regm(:, :) ! nucleation regime temperature threshold [C] + real(kind_phys), intent(out) :: subgrid_diag(:, :) ! ice nucleation subgrid saturation factor [1] + real(kind_phys), intent(out) :: trop_pd(:, :) ! chemical tropopause probability [1] + + ! Pre-existing ice diagnostics (output regardless, zeroed if not use_preexisting_ice_) + real(kind_phys), intent(out) :: fhom(:, :) ! fraction of cirrus with hom. freezing [fraction] + real(kind_phys), intent(out) :: wice(:, :) ! vertical velocity reduction from preexisting ice [m/s] + real(kind_phys), intent(out) :: weff(:, :) ! effective vertical velocity for ice nucleation [m/s] + + ! Pre-existing ice aerosol diagnostics (zeroed if not use_preexisting_ice_) + real(kind_phys), intent(out) :: INnso4(:, :) ! so4 number conc. tendency to ice nucleation [m-3 s-1] + real(kind_phys), intent(out) :: INnbc(:, :) ! bc number conc. tendency to ice nucleation [m-3 s-1] + real(kind_phys), intent(out) :: INndust(:, :) ! dust number conc. tendency to ice nucleation [m-3 s-1] + real(kind_phys), intent(out) :: INondust(:, :) ! dust number conc. tendency from ice nucleation [m-3 s-1] + real(kind_phys), intent(out) :: INhet(:, :) ! heterogeneous IN number tendency [m-3 s-1] + real(kind_phys), intent(out) :: INhom(:, :) ! homogeneous IN number tendency [m-3 s-1] + real(kind_phys), intent(out) :: INFrehom(:, :) ! homogeneous ice nucleation frequency [1] + real(kind_phys), intent(out) :: INFreIN(:, :) ! ice nucleation frequency [1] + + ! Constituent tendencies: CCPP framework applies these via tendency updater. + ! For MAM/CARMA + use_preexisting_ice: interstitial-to-cloud-borne transfer. + ! For BAM and aquaplanet: zeroed (no constituent tendencies). + real(kind_phys), intent(out) :: ptend_q(:, :, :) ! constituent tendencies [kg kg-1 s-1] (ncol, pver, pcnst) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + class(aerosol_properties), pointer :: aeroprops + class(aerosol_state), pointer :: aerostate + + integer :: i, k, m, l + integer :: idxtmp + + real(kind_phys) :: rho(ncol,pver) + real(kind_phys) :: qs(ncol) + real(kind_phys) :: es(ncol) + real(kind_phys) :: gammas(ncol) + real(kind_phys) :: relhum(ncol,pver) + real(kind_phys) :: icldm(ncol,pver) + real(kind_phys) :: subgrid(ncol,pver) + + real(kind_phys) :: dust_num_col(ncol,pver) + real(kind_phys) :: sulf_num_col(ncol,pver) + real(kind_phys) :: soot_num_col(ncol,pver) + real(kind_phys) :: sulf_num_tot_col(ncol,pver) + + real(kind_phys) :: so4_num, dst_num, soot_num + real(kind_phys) :: oso4_num, odst_num, osoot_num + real(kind_phys) :: so4_num_st_cr_tot + real(kind_phys) :: dso4_num, ramp + + ! Aerosol transfer working arrays (allocated only for MAM/CARMA) + real(kind_phys), allocatable :: size_wght(:,:,:,:) ! (ncol,pver,nbins,nmaxspc) + real(kind_phys), allocatable :: amb_num_bins(:,:,:) ! (ncol,pver,nbins) + real(kind_phys), pointer :: num_col(:,:) + real(kind_phys), pointer :: amb_mmr(:,:) + real(kind_phys), pointer :: cld_mmr(:,:) + + character(len=32) :: spectype + real(kind_phys) :: wght + real(kind_phys) :: delmmr, delmmr_sum + real(kind_phys) :: delnum, delnum_sum + + real(kind_phys), parameter :: per_cm3 = 1.e-6_kind_phys ! m-3 to cm-3 + + !------------------------------------------------------------------------------- + + errmsg = '' + errflg = 0 + + nullify(aeroprops) + nullify(aerostate) + + ! Initialize all outputs + naai(:,:) = 0._kind_phys + naai_hom(:,:) = 0._kind_phys + nihf(:,:) = 0._kind_phys + niimm(:,:) = 0._kind_phys + nidep(:,:) = 0._kind_phys + nimey(:,:) = 0._kind_phys + regm(:,:) = 0._kind_phys + trop_pd(:,:) = 0._kind_phys + fhom(:,:) = 0._kind_phys + wice(:,:) = 0._kind_phys + weff(:,:) = 0._kind_phys + ptend_q(:,:,:) = 0._kind_phys + + INnso4(:,:) = 0.0_kind_phys + INnbc(:,:) = 0.0_kind_phys + INndust(:,:) = 0.0_kind_phys + INondust(:,:) = 0.0_kind_phys + INhet(:,:) = 0.0_kind_phys + INhom(:,:) = 0.0_kind_phys + INFrehom(:,:) = 0.0_kind_phys + INFreIN(:,:) = 0.0_kind_phys + + !--------------------------------------------------------------------------- + ! Prepare abstract aerosol properties and state object. + !--------------------------------------------------------------------------- + if (iaermod_selected_ > 0) then + aeroprops => aerosol_instances_get_props(iaermod_selected_, list_idx=0) + aerostate => aerosol_instances_get_state(iaermod_selected_, list_idx=0) + end if + + ! Allocate aerosol transfer arrays if needed + if (associated(aeroprops) .and. nbins_ > 0) then + allocate(size_wght(ncol, pver, nbins_, nmaxspc_)) + allocate(amb_num_bins(ncol, pver, nbins_)) + end if + + !--------------------------------------------------------------------------- + ! Prepare derived input quantities + !--------------------------------------------------------------------------- + rho(:ncol,:) = pmid(:ncol,:) / (rair * t(:ncol,:)) + + ! Compute subgrid saturation factor from tropopause level + subgrid(:,:) = 0._kind_phys + do k = top_lev, pver + do i = 1, ncol + trop_pd(i, tropLev_chem(i)) = 1._kind_phys + + if (k <= tropLev_chem(i)) then + if (nucleate_ice_subgrid_strat_ == -1._kind_phys) then + subgrid(i, k) = 1._kind_phys / qsatfac(i, k) + else + subgrid(i, k) = nucleate_ice_subgrid_strat_ + end if + else + if (nucleate_ice_subgrid_ == -1._kind_phys) then + subgrid(i, k) = 1._kind_phys / qsatfac(i, k) + else + subgrid(i, k) = nucleate_ice_subgrid_ + end if + end if + end do + end do + subgrid_diag(:ncol,:pver) = subgrid(:ncol,:pver) + + ! Compute relative humidity and ice cloud fraction + do k = top_lev, pver + call qsat_water(t(1:ncol,k), pmid(1:ncol,k), & + es(1:ncol), qs(1:ncol), ncol, gam=gammas(1:ncol)) + do i = 1, ncol + relhum(i,k) = qv(i,k) / qs(i) + icldm(i,k) = max(aist(i,k), mincld_) + end do + end do + + ! Collect number densities [# cm-3] for dust, sulfate, and soot. + ! Unified for all aerosol models via the abstract interface. + dust_num_col = 0._kind_phys + sulf_num_col = 0._kind_phys + sulf_num_tot_col = 0._kind_phys + soot_num_col = 0._kind_phys + + if (associated(aeroprops) .and. associated(aerostate)) then + call aerostate%nuclice_get_numdens( aeroprops, use_preexisting_ice_, & + ncol, pver, rho, & + dust_num_col, sulf_num_col, & + soot_num_col, sulf_num_tot_col ) + + ! Prepare per-bin ambient numbers and size weights for aerosol transfer + do m = 1, nbins_ + call aerostate%get_ambient_num(m, num_col) + amb_num_bins(:ncol,:,m) = num_col(:ncol,:) + do l = 1, aeroprops%nspecies(m) + call aeroprops%species_type(m, l, spectype) + call aerostate%icenuc_size_wght(m, ncol, pver, spectype, & + use_preexisting_ice_, size_wght(:,:,m,l)) + end do + end do + end if + + !--------------------------------------------------------------------------- + ! Column loop: call nucleati kernel + !--------------------------------------------------------------------------- + do k = top_lev, pver + do i = 1, ncol + so4_num_st_cr_tot = 0._kind_phys + if (t(i,k) < tmelt - 5._kind_phys) then + + ! Set aerosol number for so4, soot, and dust [#/cm3] + so4_num = sulf_num_col(i,k) + dst_num = dust_num_col(i,k) + so4_num_st_cr_tot = sulf_num_tot_col(i,k) + + ! *** Turn off soot nucleation *** + ! this mod is reproduced verbatim from CAM as of cam6_4_162 + ! (hplin, 4/7/26) + soot_num = 0.0_kind_phys + + if (use_nucleati_tendencies_) then + ! PUMAS v1.21+ path (CAM7+) + call nucleati( & + wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), & + qc(i,k), qi(i,k), ni(i,k), rho(i,k), & + so4_num, dst_num, soot_num, subgrid(i,k), & + naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), & + wice(i,k), weff(i,k), fhom(i,k), regm(i,k), & + oso4_num, odst_num, osoot_num, & + call_frm_zm_in=.false., add_preexisting_ice_in=.false.) + else + ! Legacy path + ! (no optional args - defaults to: call_frm_zm = .false. + ! and add_preexisting_ice = .true) + ! this path is unsupported in CAM-SIMA (see below) (hplin, 4/7/26) + call nucleati( & + wsubi(i,k), t(i,k), pmid(i,k), relhum(i,k), icldm(i,k), & + qc(i,k), qi(i,k), ni(i,k), rho(i,k), & + so4_num, dst_num, soot_num, subgrid(i,k), & + naai(i,k), nihf(i,k), niimm(i,k), nidep(i,k), nimey(i,k), & + wice(i,k), weff(i,k), fhom(i,k), regm(i,k), & + oso4_num, odst_num, osoot_num) + end if + + !--------------------------------------------------------------------- + ! MAM / CARMA only: + ! Move aerosol used for nucleation from interstitial to cloud-borne, + ! otherwise the same coarse mode aerosols will be available again in + ! the next timestep and will suppress homogeneous freezing. + ! Replaces nucleate_ice_cam L620-692. + !--------------------------------------------------------------------- + if (clim_modal_carma_ .and. use_preexisting_ice_) then + + do m = 1, nbins_ + + if (aeroprops%icenuc_updates_num(m)) then + + if (amb_num_bins(i,k,m) > 0._kind_phys) then + delmmr_sum = 0._kind_phys + delnum_sum = 0._kind_phys + + ! iterate over the species within the bin + do l = 1, aeroprops%nspecies(m) + if (aeroprops%icenuc_updates_mmr(m, l)) then + + call aeroprops%species_type(m, l, spectype) + + wght = size_wght(i,k,m,l) + + if (wght > 0._kind_phys) then + + idxtmp = aer_cnst_idx_(m, l) + + call aerostate%get_ambient_mmr(species_ndx=l, bin_ndx=m, mmr=amb_mmr) + call aerostate%get_cldbrne_mmr(species_ndx=l, bin_ndx=m, mmr=cld_mmr) + + ! determine change in aerosol mass + delmmr = 0._kind_phys + delnum = 0._kind_phys + if (trim(spectype) == 'dust') then + if (dst_num > 0._kind_phys) then + delmmr = (odst_num / dst_num) * icldm(i,k) * amb_mmr(i,k) * wght + delnum = (odst_num * icldm(i,k)) / rho(i,k) / per_cm3 + end if + else if (trim(spectype) == 'sulfate') then + if (so4_num > 0._kind_phys) then + delmmr = (oso4_num / so4_num) * icldm(i,k) * amb_mmr(i,k) * wght + delnum = (oso4_num * icldm(i,k)) / rho(i,k) / per_cm3 + end if + end if + + if (idxtmp > 0) then + ! constituent tendency for transported species + ptend_q(i,k,idxtmp) = -delmmr / dtime + else + ! apply change of mass to not-transported species directly + amb_mmr(i,k) = amb_mmr(i,k) - delmmr + end if + cld_mmr(i,k) = cld_mmr(i,k) + delmmr + + delmmr_sum = delmmr_sum + delmmr + delnum_sum = delnum_sum + delnum + end if + end if + end do + + idxtmp = aer_cnst_idx_(m, 0) + + ! update aerosol state bin and tendency for grid box i,k + call aerostate%update_bin(m, i, k, delmmr_sum, delnum_sum, & + idxtmp, dtime, ptend_q) + + end if + + end if + end do + + end if ! clim_modal_carma_ .and. use_preexisting_ice_ + + !--------------------------------------------------------------------- + ! Polar strat SO4 enhancement + !--------------------------------------------------------------------- + ! Liu&Penner does not generate enough nucleation in the polar winter + ! stratosphere, which affects surface area density, dehydration and + ! ozone chemistry. Part of this is that there are a larger number of + ! particles in the accumulation mode than in the Aitken mode. In volcanic + ! periods, the coarse mode may also be important. As a short + ! term work around, include the accumulation and coarse mode particles + ! and assume a larger fraction of the sulfates nucleate in the polar + ! stratosphere. + ! + ! Do not include the tropopause level, as stratospheric aerosols + ! only exist above the tropopause level. + ! + ! NOTE: This may still not represent the proper particles that + ! participate in nucleation, because it doesn't include STS and NAT + ! particles. It may not represent the proper saturation threshold for + ! nucleation, and wsubi from CLUBB is probably not representative of + ! wave driven varaibility in the polar stratosphere. + if (nucleate_ice_use_troplev_ .and. clim_modal_carma_) then + if ((k < tropLev_chem(i)) .and. & + (nucleate_ice_strat_ > 0._kind_phys) .and. & + (oso4_num > 0._kind_phys)) then + dso4_num = max(0._kind_phys, & + (nucleate_ice_strat_ * so4_num_st_cr_tot - oso4_num) & + * 1e6_kind_phys / rho(i,k)) + naai(i,k) = naai(i,k) + dso4_num + nihf(i,k) = nihf(i,k) + dso4_num + end if + else + ! Fallback: pressure/latitude heuristic (maintains backwards compatibility) + if (pmid(i,k) <= 12500._kind_phys .and. & + pmid(i,k) > 100._kind_phys .and. & + abs(lat(i)) >= 60._kind_phys * pi / 180._kind_phys) then + ramp = 1._kind_phys - min(1._kind_phys, & + max(0._kind_phys, (pmid(i,k) - 10000._kind_phys) / 2500._kind_phys)) + + if (oso4_num > 0._kind_phys) then + dso4_num = (max(oso4_num, ramp * nucleate_ice_strat_ * so4_num) & + - oso4_num) * 1e6_kind_phys / rho(i,k) + naai(i,k) = naai(i,k) + dso4_num + nihf(i,k) = nihf(i,k) + dso4_num + end if + end if + end if + + !--------------------------------------------------------------------- + ! Convert to tendencies for PUMAS (CAM7+) + !--------------------------------------------------------------------- + if (use_nucleati_tendencies_) then + ! ^^ this flag was previously CAM7+ + ! PUMAS v1.21+ path: convert to tendencies [1/kg/s] and [1/m3/s] + + naai_hom(i,k) = nihf(i,k) / dtime + naai(i,k) = naai(i,k) / dtime + + ! output activated ice (convert from # kg-1 to # m-3 s-1) + nihf(i,k) = nihf(i,k) * rho(i,k) / dtime + niimm(i,k) = niimm(i,k) * rho(i,k) / dtime + nidep(i,k) = nidep(i,k) * rho(i,k) / dtime + nimey(i,k) = nimey(i,k) * rho(i,k) / dtime + + if (use_preexisting_ice_) then + INnso4(i,k) = so4_num * 1e6_kind_phys / dtime ! # cm-3 -> # m-3 s-1 + INnbc(i,k) = soot_num * 1e6_kind_phys / dtime + INndust(i,k) = dst_num * 1e6_kind_phys / dtime + INondust(i,k)= odst_num * 1e6_kind_phys / dtime + INFreIN(i,k) = 1.0_kind_phys ! 1, ice nucleation occurred + INhet(i,k) = (niimm(i,k) + nidep(i,k)) ! # m-3 s-1, nimey not in cirrus + INhom(i,k) = nihf(i,k) ! # m-3 s-1 + if (INhom(i,k) > 1e3_kind_phys) then ! > 1/L + INFrehom(i,k) = 1.0_kind_phys ! 1, hom freezing occurred + end if + + ! exclude no ice nucleation + if ((INFrehom(i,k) < 0.5_kind_phys) .and. & + (INhet(i,k) < 1.0_kind_phys)) then + INnso4(i,k) = 0.0_kind_phys + INnbc(i,k) = 0.0_kind_phys + INndust(i,k) = 0.0_kind_phys + INondust(i,k) = 0.0_kind_phys + INFreIN(i,k) = 0.0_kind_phys + INhet(i,k) = 0.0_kind_phys + INhom(i,k) = 0.0_kind_phys + INFrehom(i,k) = 0.0_kind_phys + wice(i,k) = 0.0_kind_phys + weff(i,k) = 0.0_kind_phys + fhom(i,k) = 0.0_kind_phys + end if + end if + else + ! Legacy path (CAM5/6): unsupported in CAM-SIMA CCPP. + ! Kept commented for reference. If needed, uncomment and + ! change output units from tendencies to concentrations. + ! + ! naai_hom(i,k) = nihf(i,k) + ! nihf(i,k) = nihf(i,k) * rho(i,k) + ! niimm(i,k) = niimm(i,k) * rho(i,k) + ! nidep(i,k) = nidep(i,k) * rho(i,k) + ! nimey(i,k) = nimey(i,k) * rho(i,k) + + errflg = 1 + errmsg = 'non-CAM7 path for nucleate_ice is unsupported in CCPP' + return + end if ! use_nucleati_tendencies_ + end if ! T < tmelt - 5 + end do ! i + end do ! k + + ! Clean up + if (allocated(size_wght)) deallocate(size_wght) + if (allocated(amb_num_bins)) deallocate(amb_num_bins) + + end subroutine nucleate_ice_ccpp_run + +end module nucleate_ice_ccpp diff --git a/schemes/microp_aero/ccpp/nucleate_ice_ccpp.meta b/schemes/microp_aero/ccpp/nucleate_ice_ccpp.meta new file mode 100644 index 00000000..93995c63 --- /dev/null +++ b/schemes/microp_aero/ccpp/nucleate_ice_ccpp.meta @@ -0,0 +1,350 @@ +[ccpp-table-properties] + name = nucleate_ice_ccpp + type = scheme + dependencies = ../nucleate_ice.F90 + +[ccpp-arg-table] + name = nucleate_ice_ccpp_init + type = scheme +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ const_props ] + standard_name = ccpp_constituent_properties + units = none + type = ccpp_constituent_prop_ptr_t + dimensions = (number_of_ccpp_constituents) + intent = in +[ use_preexisting_ice ] + standard_name = do_preexisting_ice_treatment_in_ice_nucleation + units = flag + type = logical + dimensions = () + intent = in +[ use_hetfrz_classnuc ] + standard_name = do_heterogeneous_freezing_by_classical_nucleation_theory + units = flag + type = logical + dimensions = () + intent = in +[ nucleate_ice_incloud ] + standard_name = do_incloud_ice_nucleation + units = flag + type = logical + dimensions = () + intent = in +[ pi ] + standard_name = pi_constant + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ nucleate_ice_subgrid ] + standard_name = tunable_parameter_for_subgrid_saturation_scaling_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ nucleate_ice_subgrid_strat ] + standard_name = tunable_parameter_for_subgrid_saturation_scaling_in_stratosphere_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ nucleate_ice_strat ] + standard_name = tunable_parameter_for_polar_stratospheric_sulfate_nucleation_fraction_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ nucleate_ice_use_troplev ] + standard_name = do_tropopause_level_based_stratospheric_adjustment_for_ice_nucleation + units = flag + type = logical + dimensions = () + intent = in +[ prog_modal_aero ] + standard_name = do_prognostic_modal_aerosol + units = flag + type = logical + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = nucleate_ice_ccpp_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pcnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qv ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ qc ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ qi ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ ni ] + standard_name = mass_number_concentration_of_ice_wrt_moist_air_and_condensed_water + units = kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = True +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ pi ] + standard_name = pi_constant + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ wsubi ] + standard_name = subgrid_scale_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ aist ] + standard_name = stratiform_cloud_ice_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ tropLev_chem ] + standard_name = tropopause_vertical_layer_index_from_chemical_method + units = index + type = integer + dimensions = (horizontal_loop_extent) + intent = in +[ qsatfac ] + # pbuf qsatfac + standard_name = subgrid_cloud_water_saturation_scaling_factor_for_microphysics + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ lat ] + standard_name = latitude + units = rad + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ naai ] + # these are tendencies for PUMAS v1.21+ (CAM7+) which I assume is the only supported path for SIMA. + standard_name = tendency_of_number_concentration_of_activated_ice_nuclei + units = kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ naai_hom ] + # these are tendencies for PUMAS v1.21+ (CAM7+) which I assume is the only supported path for SIMA. + standard_name = tendency_of_number_concentration_of_activated_ice_nuclei_due_to_homogenous_freezing + units = kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ nihf ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_homogeneous_freezing + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ niimm ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_immersion_freezing + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ nidep ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_deposition_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ nimey ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_meyers_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ regm ] + standard_name = ice_nucleation_regime_temperature_threshold + units = degC + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ subgrid_diag ] + standard_name = subgrid_saturation_scaling_factor_diagnostic_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ trop_pd ] + # in A flag is accummulated to be a probability distribution + standard_name = vertical_layer_indicator_of_chemical_tropopause + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fhom ] + standard_name = fraction_of_cirrus_with_homogeneous_freezing + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ wice ] + standard_name = vertical_velocity_reduction_from_preexisting_ice_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ weff ] + standard_name = effective_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INnso4 ] + standard_name = tendency_of_sulfate_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INnbc ] + standard_name = tendency_of_black_carbon_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INndust ] + standard_name = tendency_of_dust_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INondust ] + standard_name = tendency_of_dust_number_concentration_output_from_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INhet ] + standard_name = tendency_of_heterogeneous_ice_nuclei_number_concentration + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INhom ] + standard_name = tendency_of_homogeneous_ice_nuclei_number_concentration + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INFrehom ] + standard_name = frequency_of_homogeneous_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ INFreIN ] + standard_name = frequency_of_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ ptend_q ] + standard_name = ccpp_constituent_tendencies + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/ccpp/nucleate_ice_ccpp_namelist.xml b/schemes/microp_aero/ccpp/nucleate_ice_ccpp_namelist.xml new file mode 100644 index 00000000..7be2c85a --- /dev/null +++ b/schemes/microp_aero/ccpp/nucleate_ice_ccpp_namelist.xml @@ -0,0 +1,114 @@ + + + + + + + + logical + microphys + nucleate_ice_nl + do_preexisting_ice_treatment_in_ice_nucleation + flag + + Switch to turn on treatment of pre-existing ice in the ice nucleation code. + Default: TRUE (CAM6 and CAM7) + + + .true. + .false. + + + + + logical + microphys + nucleate_ice_nl + do_incloud_ice_nucleation + flag + + Switch to determine whether ice nucleation happens using the incloud (true) or + the gridbox average (false) relative humidity. When true, it is assumed that + the incloud relative humidity for nucleation is 1. + Default: false. + + + .false. + .true. + + + + + real + microphys + nucleate_ice_nl + tunable_parameter_for_subgrid_saturation_scaling_for_ice_nucleation + 1 + + Subgrid scaling factor for relative humidity in ice nucleation code. If it has + a value of -1, then indicates that the subgrid scaling factor will be + calculated on the fly as 1 / qsatfac (i.e. the saturation scaling factor). + Default: 1.2 (MG1/MG2/MG3) + + + 1.2D0 + + + + + real + microphys + nucleate_ice_nl + tunable_parameter_for_subgrid_saturation_scaling_in_stratosphere_for_ice_nucleation + 1 + + Subgrid scaling factor for relative humidity in ice nucleation code in the + stratosphere. If it has a value of -1, then indicates that the subgrid + scaling factor will be calculated on the fly as 1 / qsatfac (i.e. the + saturation scaling factor). + Default: 1.2 (MG1/MG2/MG3) + + + 1.2D0 + + + + + real + microphys + nucleate_ice_nl + tunable_parameter_for_polar_stratospheric_sulfate_nucleation_fraction_for_ice_nucleation + 1 + + Fraction of Aitken mode sulfate particles assumed to nucleate ice in the polar + stratosphere. Provides an increase in homogeneous freezing over the Liu and Penner + method. + + Temporary solution to adjust ice surface area density and dehydration in + the polar stratosphere where there doesn't seem to be enough nucleation. A value + of 1.0 disables the effect. + Default: 1.0 + + + 1.0D0 + + + + + logical + microphys + nucleate_ice_nl + do_tropopause_level_based_stratospheric_adjustment_for_ice_nucleation + flag + + Indicates whether to use the tropopause level to determine where to adjust + nucleation for the stratosphere (true) or whether to use a hard coded transition + level from 100 to 125 hPa applied only in the polar regions (false). + Default: true + + + .true. + + + + diff --git a/schemes/microp_aero/compute_subgrid_vertical_velocity.F90 b/schemes/microp_aero/compute_subgrid_vertical_velocity.F90 new file mode 100644 index 00000000..a2b35bf9 --- /dev/null +++ b/schemes/microp_aero/compute_subgrid_vertical_velocity.F90 @@ -0,0 +1,145 @@ +! Compute raw subgrid-scale vertical velocity from TKE, KVH, or CLUBB WP2. +module compute_subgrid_vertical_velocity + + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: compute_subgrid_vertical_velocity_tke_run + public :: compute_subgrid_vertical_velocity_kvh_run + public :: compute_subgrid_vertical_velocity_clubb_run + +contains + +!> \section arg_table_compute_subgrid_vertical_velocity_tke_run Argument Table +!! \htmlinclude compute_subgrid_vertical_velocity_tke_run.html + + ! Note: despite the tke name, the "tke" used here is the CAM5 diag_TKE scheme TKE + ! and is not the TKE originating from CLUBB, and thus this scheme is not used for + ! CAM6+, it is used for CAM5 only. + ! For CAM6+ do not use TKE, derive from WP2_nadv - see the clubb variant. hplin 4/20/26 + subroutine compute_subgrid_vertical_velocity_tke_run( & + ncol, pver, top_lev, & + tke, & + wsub, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev ! top vertical level for cloud physics + + real(kind_phys), intent(in) :: tke(:, :) ! turbulent kinetic energy at interfaces [m2 s-2] + + ! Output arguments + real(kind_phys), intent(out) :: wsub(:, :) ! subgrid vertical velocity [m s-1] + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + + errmsg = '' + errflg = 0 + + wsub(:, :top_lev-1) = 0._kind_phys + + do k = top_lev, pver + do i = 1, ncol + wsub(i,k) = sqrt(0.5_kind_phys * (tke(i,k) + tke(i,k+1)) * (2._kind_phys / 3._kind_phys)) + wsub(i,k) = min(wsub(i,k), 10._kind_phys) + end do + end do + + end subroutine compute_subgrid_vertical_velocity_tke_run + +!> \section arg_table_compute_subgrid_vertical_velocity_kvh_run Argument Table +!! \htmlinclude compute_subgrid_vertical_velocity_kvh_run.html + subroutine compute_subgrid_vertical_velocity_kvh_run( & + ncol, pver, top_lev, & + kvh, & + wsub, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev ! top vertical level for cloud physics + + real(kind_phys), intent(in) :: kvh(:, :) ! eddy diffusivity for heat at interfaces [m2 s-1] + + ! Output arguments + real(kind_phys), intent(out) :: wsub(:, :) ! subgrid vertical velocity [m s-1] + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + + errmsg = '' + errflg = 0 + + wsub(:, :top_lev-1) = 0._kind_phys + + ! Get sub-grid vertical velocity from diffusion coefficient. + ! Following Morrison et al. 2005, JAS. + ! Assume mixing length of 30 m. + ! Use maximum sub-grid vertical vel of 10 m/s. + do k = top_lev, pver + do i = 1, ncol + wsub(i,k) = (kvh(i,k) + kvh(i,k+1)) / 2._kind_phys / 30._kind_phys + wsub(i,k) = min(wsub(i,k), 10._kind_phys) + end do + end do + + end subroutine compute_subgrid_vertical_velocity_kvh_run + +!> \section arg_table_compute_subgrid_vertical_velocity_clubb_run Argument Table +!! \htmlinclude compute_subgrid_vertical_velocity_clubb_run.html + subroutine compute_subgrid_vertical_velocity_clubb_run( & + ncol, pver, pverp, top_lev, & + wp2, & + wsub, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + integer, intent(in) :: top_lev ! top vertical level for cloud physics + + real(kind_phys), intent(in) :: wp2(:, :) ! CLUBB variance of vertical velocity at interfaces [m2 s-2] + + ! Output arguments + real(kind_phys), intent(out) :: wsub(:, :) ! subgrid vertical velocity [m s-1] + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + real(kind_phys) :: tke(ncol, pver+1) + + errmsg = '' + errflg = 0 + + ! Convert wp2 to TKE: tke = (3/2) * wp2 from [pver+1, top_lev] + ! This matches CAM microp_aero.F90 CLUBB_SGS branch exactly. + tke(:ncol,top_lev:pverp) = (3._kind_phys/2._kind_phys)*wp2(:ncol,1:pverp-top_lev+1) + tke(:ncol,1:top_lev-1) = 0._kind_phys + + wsub(:, :top_lev-1) = 0._kind_phys + + do k = top_lev, pver + do i = 1, ncol + wsub(i,k) = sqrt(0.5_kind_phys * (tke(i,k) + tke(i,k+1)) * (2._kind_phys / 3._kind_phys)) + wsub(i,k) = min(wsub(i,k), 10._kind_phys) + end do + end do + + end subroutine compute_subgrid_vertical_velocity_clubb_run + +end module compute_subgrid_vertical_velocity diff --git a/schemes/microp_aero/compute_subgrid_vertical_velocity.meta b/schemes/microp_aero/compute_subgrid_vertical_velocity.meta new file mode 100644 index 00000000..aade6c8e --- /dev/null +++ b/schemes/microp_aero/compute_subgrid_vertical_velocity.meta @@ -0,0 +1,155 @@ +[ccpp-table-properties] + name = compute_subgrid_vertical_velocity_tke + type = scheme + +[ccpp-arg-table] + name = compute_subgrid_vertical_velocity_tke_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ tke ] + standard_name = specific_turbulent_kinetic_energy_at_interfaces + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = compute_subgrid_vertical_velocity_clubb + type = scheme + +[ccpp-arg-table] + name = compute_subgrid_vertical_velocity_clubb_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pverp ] + standard_name = vertical_interface_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ wp2 ] + standard_name = advected_variance_of_vertical_velocity_at_interfaces + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = compute_subgrid_vertical_velocity_kvh + type = scheme + +[ccpp-arg-table] + name = compute_subgrid_vertical_velocity_kvh_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ kvh ] + standard_name = eddy_heat_diffusivity_at_interfaces + units = m2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/dust_default_radii.F90 b/schemes/microp_aero/dust_default_radii.F90 new file mode 100644 index 00000000..4dd105bb --- /dev/null +++ b/schemes/microp_aero/dust_default_radii.F90 @@ -0,0 +1,40 @@ +! Set default dust radii for 4 size bins used in contact freezing. +! For modal aerosols, bin 3 is overwritten later with wet diameter. +module dust_default_radii + + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: dust_default_radii_timestep_init + + real(kind_phys), parameter :: rn_dst1 = 0.258e-6_kind_phys + real(kind_phys), parameter :: rn_dst2 = 0.717e-6_kind_phys + real(kind_phys), parameter :: rn_dst3 = 1.576e-6_kind_phys + real(kind_phys), parameter :: rn_dst4 = 3.026e-6_kind_phys + +contains + +!> \section arg_table_dust_default_radii_timestep_init Argument Table +!! \htmlinclude dust_default_radii_timestep_init.html + subroutine dust_default_radii_timestep_init(ncol, pver, ndust, rndst, errmsg, errflg) + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: ndust + + real(kind_phys), intent(out) :: rndst(:,:,:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + rndst(:ncol, :pver, 1) = rn_dst1 + rndst(:ncol, :pver, 2) = rn_dst2 + rndst(:ncol, :pver, 3) = rn_dst3 + rndst(:ncol, :pver, 4) = rn_dst4 + + end subroutine dust_default_radii_timestep_init + +end module dust_default_radii diff --git a/schemes/microp_aero/dust_default_radii.meta b/schemes/microp_aero/dust_default_radii.meta new file mode 100644 index 00000000..785d4394 --- /dev/null +++ b/schemes/microp_aero/dust_default_radii.meta @@ -0,0 +1,43 @@ +[ccpp-table-properties] + name = dust_default_radii + type = scheme + +[ccpp-arg-table] + name = dust_default_radii_timestep_init + type = scheme +[ ncol ] + standard_name = horizontal_dimension + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ ndust ] + standard_name = dust_size_bin_dimension + units = count + type = integer + dimensions = () + intent = in +[ rndst ] + standard_name = dust_radii_by_size_bin + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension, dust_size_bin_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/ndrop_bam.F90 b/schemes/microp_aero/ndrop_bam.F90 new file mode 100644 index 00000000..bc4ad403 --- /dev/null +++ b/schemes/microp_aero/ndrop_bam.F90 @@ -0,0 +1,628 @@ +! Droplet activation by bulk aerosols. +module ndrop_bam + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: ndrop_bam_init + public :: ndrop_bam_calc + + ! these are currently public parameters for use by CAM. + public :: aername, psat, ccn_name + + ! module parameters: + + ! # of supersaturations to calc ccn concentration + integer, parameter :: psat = 6 + ! supersaturation (%) to determine ccn concentration + real(kind_phys), parameter :: supersat(psat) = & + (/0.02_kind_phys, 0.05_kind_phys, 0.1_kind_phys, 0.2_kind_phys, 0.5_kind_phys, 1.0_kind_phys/) + ! corresponding diagnostic names: + character(len=4), parameter :: ccn_name(psat) = & + (/'CCN1', 'CCN2', 'CCN3', 'CCN4', 'CCN5', 'CCN6'/) + real(kind_phys) :: super(psat) + + real(kind_phys), allocatable :: ccnfact(:, :) + real(kind_phys), allocatable :: alogsig(:) ! natl log of geometric standard dev of aerosol + real(kind_phys), allocatable :: exp45logsig(:) + real(kind_phys), allocatable :: argfactor(:) + real(kind_phys), allocatable :: amcube(:) ! cube of dry mode radius (m) + real(kind_phys), allocatable :: smcrit(:) ! critical supersaturation for activation + real(kind_phys), allocatable :: lnsm(:) ! ln(smcrit) + real(kind_phys), allocatable :: amcubefactor(:) ! factors for calculating mode radius + real(kind_phys), allocatable :: smcritfactor(:) ! factors for calculating critical supersaturation + real(kind_phys), allocatable :: f1(:), f2(:) ! abdul-razzak functions of width + + real(kind_phys) :: pi ! pi + real(kind_phys) :: aten + real(kind_phys) :: third, sixth + real(kind_phys) :: sq2, sqpi + real(kind_phys) :: alogten, alog2, alog3, alogaten + + ! aerosol properties + character(len=20), allocatable :: aername(:) + real(kind_phys), allocatable :: dryrad_aer(:) + real(kind_phys), allocatable :: density_aer(:) + real(kind_phys), allocatable :: hygro_aer(:) + real(kind_phys), allocatable :: dispersion_aer(:) + real(kind_phys), allocatable :: num_to_mass_aer(:) + + integer :: naer_all ! number of aerosols affecting climate + integer :: idxsul = -1 ! index in aerosol list for sulfate + integer :: idxdst2 = -1 ! index in aerosol list for dust2 + integer :: idxdst3 = -1 ! index in aerosol list for dust3 + integer :: idxdst4 = -1 ! index in aerosol list for dust4 + +contains + + ! Initialize constants for droplet activation by bulk aerosols + subroutine ndrop_bam_init(amIRoot, iulog, mwh2o, r_universal, tmelt, rhoh2o, & + naer_all_out) + + use shr_spfn_mod, only: erf => shr_spfn_erf + + use aerosol_instances_mod, only: aerosol_instances_get_props, aerosol_instances_get_num_models + use aerosol_properties_mod, only: aerosol_properties + + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + real(kind_phys), intent(in) :: mwh2o ! molecular weight of water (kg/kmol) + real(kind_phys), intent(in) :: r_universal ! universal gas constant (J/K/kmol) + real(kind_phys), intent(in) :: tmelt ! freezing point of water (K) + real(kind_phys), intent(in) :: rhoh2o ! density of liquid water (kg/m3) + integer, intent(out) :: naer_all_out + + integer :: l, m, iaer, iaermod + real(kind_phys) :: surften ! surface tension of water w/respect to air (N/m) + real(kind_phys) :: arg + class(aerosol_properties), pointer :: aero_props_bam + !------------------------------------------------------------------------------- + + ! Access the physical properties of the bulk aerosols that are affecting the climate + ! by using the abstract aerosol properties interface. + + ! Find BAM properties object from factory + aero_props_bam => null() + do iaermod = 1, aerosol_instances_get_num_models() + aero_props_bam => aerosol_instances_get_props(iaermod, 0) + if (associated(aero_props_bam)) then + if (aero_props_bam%model_is('BAM')) exit + end if + aero_props_bam => null() + end do + + if (associated(aero_props_bam)) then + naer_all = aero_props_bam%nbins() + else + naer_all = 0 + end if + + allocate ( & + aername(naer_all), & + dryrad_aer(naer_all), & + density_aer(naer_all), & + hygro_aer(naer_all), & + dispersion_aer(naer_all), & + num_to_mass_aer(naer_all)) + + do iaer = 1, naer_all + call aero_props_bam%get(iaer, 1, & + specname=aername(iaer), & + dryrad=dryrad_aer(iaer), & + density=density_aer(iaer), & + hygro=hygro_aer(iaer), & + num_to_mass_aer=num_to_mass_aer(iaer)) + dispersion_aer(iaer) = exp(aero_props_bam%alogsig(iaer)) + + ! Look for sulfate and dust aerosols in this list (Bulk aerosol only) + if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer + if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer + if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer + if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer + + end do + + if (amIRoot) then + write (iulog, *) 'ndrop_bam_init: iaer, name, dryrad, density, hygro, dispersion, num_to_mass' + do iaer = 1, naer_all + write (iulog, *) iaer, aername(iaer), dryrad_aer(iaer), density_aer(iaer), hygro_aer(iaer), & + dispersion_aer(iaer), num_to_mass_aer(iaer) + end do + if (idxsul < 1) then + write (iulog, *) 'ndrop_bam_init: SULFATE aerosol properties NOT FOUND' + else + write (iulog, *) 'ndrop_bam_init: SULFATE aerosol properties FOUND at index ', idxsul + end if + end if + + ! set parameters for droplet activation, + ! following abdul-razzak and ghan 2000, JGR + third = 1._kind_phys/3._kind_phys + sixth = 1._kind_phys/6._kind_phys + sq2 = sqrt(2._kind_phys) + pi = 4._kind_phys*atan(1.0_kind_phys) ! note: this should probably use physconst pi but here for b4b + sqpi = sqrt(pi) + surften = 0.076_kind_phys + aten = 2._kind_phys*mwh2o*surften/(r_universal*tmelt*rhoh2o) + alogaten = log(aten) + alog2 = log(2._kind_phys) + alog3 = log(3._kind_phys) + super(:) = 0.01_kind_phys*supersat(:) + + allocate ( & + alogsig(naer_all), & + exp45logsig(naer_all), & + argfactor(naer_all), & + f1(naer_all), & + f2(naer_all), & + amcubefactor(naer_all), & + smcritfactor(naer_all), & + amcube(naer_all), & + smcrit(naer_all), & + lnsm(naer_all), & + ccnfact(psat, naer_all)) + + do m = 1, naer_all + + ! Skip aerosols that don't have a dispersion defined. + if (dispersion_aer(m) == 0._kind_phys) cycle + + alogsig(m) = log(dispersion_aer(m)) + exp45logsig(m) = exp(4.5_kind_phys*alogsig(m)*alogsig(m)) + argfactor(m) = 2._kind_phys/(3._kind_phys*sqrt(2._kind_phys)*alogsig(m)) + f1(m) = 0.5_kind_phys*exp(2.5_kind_phys*alogsig(m)*alogsig(m)) + f2(m) = 1._kind_phys + 0.25_kind_phys*alogsig(m) + amcubefactor(m) = 3._kind_phys/(4._kind_phys*pi*exp45logsig(m)*density_aer(m)) + smcritfactor(m) = 2._kind_phys*aten*sqrt(aten/(27._kind_phys*max(1.e-10_kind_phys, hygro_aer(m)))) + amcube(m) = amcubefactor(m)/num_to_mass_aer(m) + + if (hygro_aer(m) .gt. 1.e-10_kind_phys) then + smcrit(m) = smcritfactor(m)/sqrt(amcube(m)) + else + smcrit(m) = 100._kind_phys + end if + lnsm(m) = log(smcrit(m)) + + do l = 1, psat + arg = argfactor(m)*log(smcrit(m)/super(l)) + if (arg < 2) then + if (arg < -2) then + ccnfact(l, m) = 1.e-6_kind_phys + else + ccnfact(l, m) = 1.e-6_kind_phys*0.5_kind_phys*erfc(arg) + end if + else + ccnfact(l, m) = 0._kind_phys + end if + end do + + end do + + naer_all_out = naer_all + + end subroutine ndrop_bam_init + + ! BAM droplet activation, contact freezing dust, and CCN diagnostics. + ! Computes naer2/maerosol from aerosol state internally via get_bulk_num_and_mass. + subroutine ndrop_bam_calc( & + aero_state, aero_props, & + ncol, pver, top_lev, & + gravit, rair, tmelt, cpair, rh2o, rhoh2o, latvap, & + rho, tair, wsub, qcld, qsmall_in, ast, numliq, deltatin, & + npccn, nacon, ccn, naer2_diag, & + errmsg, errflg) + + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state + use bulk_aerosol_state_mod, only: bulk_aerosol_state + + class(aerosol_state), intent(in) :: aero_state + class(aerosol_properties), intent(in) :: aero_props + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev + real(kind_phys), intent(in) :: gravit ! gravitational acceleration (m/s2) + real(kind_phys), intent(in) :: rair ! dry air gas constant (J/K/kg) + real(kind_phys), intent(in) :: tmelt ! freezing point of water (K) + real(kind_phys), intent(in) :: cpair ! specific heat of dry air (J/K/kg) + real(kind_phys), intent(in) :: rh2o ! water vapor gas constant (J/K/kg) + real(kind_phys), intent(in) :: rhoh2o ! density of liquid water (kg/m3) + real(kind_phys), intent(in) :: latvap ! latent heat of vaporization (J/kg) + real(kind_phys), intent(in) :: rho(:, :) ! air density (kg/m3) + real(kind_phys), intent(in) :: tair(:, :) ! temperature (K) + real(kind_phys), intent(in) :: wsub(:, :) ! sub-grid vertical velocity (m/s) + real(kind_phys), intent(in) :: qcld(:, :) ! cloud liquid mmr (kg/kg) + real(kind_phys), intent(in) :: qsmall_in ! minimum cloud liquid threshold + real(kind_phys), intent(in) :: ast(:, :) ! stratiform_cloud_area_fraction [fraction] + real(kind_phys), intent(in) :: numliq(:, :) ! droplet number (#/kg) + real(kind_phys), intent(in) :: deltatin ! timestep (s) + real(kind_phys), intent(out) :: npccn(:, :) ! droplet number tendency (#/kg/s) + real(kind_phys), intent(out) :: nacon(:, :, :) ! contact freezing dust (#/m3), (ncol,pver,4) + real(kind_phys), intent(out) :: ccn(:, :, :) ! CCN at 6 supersaturations (#/cm3), (ncol,pver,psat) + real(kind_phys), intent(out) :: naer2_diag(:, :, :) ! aerosol number conc (#/m3), (ncol,pver,naer_all) + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! minimum allowed cloud fraction + real(kind_phys), parameter :: mincld = 0.0001_kind_phys + + ! local workspace + real(kind_phys), allocatable :: naer2(:, :, :) ! aerosol number concentration [m-3] + real(kind_phys), allocatable :: maerosol(:, :, :) ! aerosol mass conc [kg m-3] + real(kind_phys) :: lcldm(ncol, pver) + real(kind_phys) :: nact + integer :: i, k, m + + errmsg = '' + errflg = 0 + + allocate (naer2(ncol, pver, naer_all), stat=errflg, errmsg=errmsg) + allocate (maerosol(ncol, pver, naer_all), stat=errflg, errmsg=errmsg) + if (errflg /= 0) return + + ! Compute aerosol number and mass concentrations from aerosol state. + ! b4b operation order: (mmr * rho) * ntm [* 2.0 for sulfate]. + select type (bam => aero_state) + type is (bulk_aerosol_state) + do m = 1, naer_all + call bam%get_bulk_num_and_mass(m, ncol, rho, naer2(:, :, m), maerosol(:, :, m)) + end do + class default + errmsg = 'ndrop_bam_calc: aero_state must be bulk_aerosol_state' + errflg = 1 + return + end select + + do k = top_lev, pver + do i = 1, ncol + lcldm(i, k) = max(ast(i, k), mincld) + end do + end do + + ! Droplet activation + npccn(:, :) = 0._kind_phys + do k = top_lev, pver + do i = 1, ncol + if (naer_all > 0 .and. qcld(i, k) >= qsmall_in) then + call activate(wsub(i, k), tair(i, k), rho(i, k), naer2(i, k, :), & + naer_all, naer_all, maerosol(i, k, :), & + gravit, rair, tmelt, cpair, rh2o, rhoh2o, latvap, & + nact, errmsg, errflg) + if (errflg /= 0) return + else + nact = 0._kind_phys + end if + npccn(i, k) = (nact*lcldm(i, k) - numliq(i, k))/deltatin + end do + end do + + ! Contact freezing: dust number concentrations for bins 2-4 + nacon(:, :, :) = 0._kind_phys + do k = top_lev, pver + do i = 1, ncol + if (tair(i, k) < 269.15_kind_phys) then + if (idxdst2 > 0) nacon(i, k, 2) = naer2(i, k, idxdst2) + if (idxdst3 > 0) nacon(i, k, 3) = naer2(i, k, idxdst3) + if (idxdst4 > 0) nacon(i, k, 4) = naer2(i, k, idxdst4) + end if + end do + end do + + ! CCN diagnostics + call ccn_diag(ncol, pver, top_lev, maerosol, naer2, ccn) + + ! Copy naer2 for diagnostic output + naer2_diag(:ncol, :, :) = naer2(:ncol, :, :) + + deallocate (naer2, maerosol) + + end subroutine ndrop_bam_calc + +!=============================================================================== + + subroutine activate( & + wbar, tair, rhoair, na, pmode, & + nmode, ma, & + gravit, rair, tmelt, cpair, rh2o, rhoh2o, latvap, & + nact, errmsg, errflg) + + use shr_spfn_mod, only: erf => shr_spfn_erf + use wv_saturation, only: qsat + + ! calculates number fraction of aerosols activated as CCN + ! assumes an internal mixture within each of up to pmode multiple aerosol modes + ! a gaussian spectrum of updrafts can be treated. + + ! mks units + + ! Abdul-Razzak and Ghan, A parameterization of aerosol activation. + ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844. + ! https://doi.org/10.1029/1999JD901161 + + ! input + integer, intent(in) :: pmode ! dimension of modes + integer, intent(in) :: nmode ! number of aerosol modes + real(kind_phys), intent(in) :: wbar ! grid cell mean vertical velocity (m/s) + real(kind_phys), intent(in) :: tair ! air temperature (K) + real(kind_phys), intent(in) :: rhoair ! air density (kg/m3) + real(kind_phys), intent(in) :: na(pmode) ! aerosol number concentration (1/m3) + real(kind_phys), intent(in) :: ma(pmode) ! aerosol mass concentration (kg/m3) + real(kind_phys), intent(in) :: gravit ! gravitational acceleration (m/s2) + real(kind_phys), intent(in) :: rair ! dry air gas constant (J/K/kg) + real(kind_phys), intent(in) :: tmelt ! freezing point of water (K) + real(kind_phys), intent(in) :: cpair ! specific heat of dry air (J/K/kg) + real(kind_phys), intent(in) :: rh2o ! water vapor gas constant (J/K/kg) + real(kind_phys), intent(in) :: rhoh2o ! density of liquid water (kg/m3) + real(kind_phys), intent(in) :: latvap ! latent heat of vaporization (J/kg) + + ! output + real(kind_phys), intent(out) :: nact ! number fraction of aerosols activated + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! local variables + integer :: maxmodes + + real(kind_phys), allocatable :: volc(:) ! total aerosol volume concentration (m3/m3) + real(kind_phys), allocatable :: eta(:) + real(kind_phys), allocatable :: smc(:) + real(kind_phys), allocatable :: etafactor2(:) + real(kind_phys), allocatable :: amcubeloc(:) + real(kind_phys), allocatable :: lnsmloc(:) + + real(kind_phys), parameter :: pref = 1013.25e2_kind_phys ! reference pressure [Pa] + + real(kind_phys) :: pres ! pressure (Pa) + real(kind_phys) :: diff0 + real(kind_phys) :: conduct0 ! thermal conductivity (Joule/m/sec/deg) + real(kind_phys) :: qs ! water vapor saturation mixing ratio + real(kind_phys) :: dqsdt ! change in qs with temperature + real(kind_phys) :: gloc ! thermodynamic function (m2/s) + real(kind_phys) :: zeta + real(kind_phys) :: lnsmax ! ln(smax) + real(kind_phys) :: alpha + real(kind_phys) :: gammaloc + real(kind_phys) :: beta + real(kind_phys) :: sqrtg + real(kind_phys) :: wnuc + real(kind_phys) :: alw + real(kind_phys) :: sqrtalw + real(kind_phys) :: smax + real(kind_phys) :: x + real(kind_phys) :: etafactor1 + real(kind_phys) :: etafactor2max + real(kind_phys) :: es + integer :: m + + errflg = 0 + errmsg = '' + + maxmodes = naer_all + allocate ( & + volc(maxmodes), & + eta(maxmodes), & + smc(maxmodes), & + etafactor2(maxmodes), & + amcubeloc(maxmodes), & + lnsmloc(maxmodes), & + stat=errflg, errmsg=errmsg) + if (errflg /= 0) return + + if (maxmodes < pmode) then + write (errmsg, *) 'ndrop_bam activate: maxmodes,pmode=', maxmodes, pmode + errflg = 1 + return + end if + + nact = 0._kind_phys + + if (nmode .eq. 1 .and. na(1) .lt. 1.e-20_kind_phys) return + + if (wbar .le. 0._kind_phys) return + + pres = rair*rhoair*tair + diff0 = 0.211e-4_kind_phys*(pref/pres)*(tair/tmelt)**1.94_kind_phys + conduct0 = (5.69_kind_phys + 0.017_kind_phys*(tair - tmelt))*4.186e2_kind_phys*1.e-5_kind_phys ! convert to J/m/s/deg + call qsat(tair, pres, es, qs) + dqsdt = latvap/(rh2o*tair*tair)*qs + alpha = gravit*(latvap/(cpair*rh2o*tair*tair) - 1._kind_phys/(rair*tair)) + gammaloc = (1 + latvap/cpair*dqsdt)/(rhoair*qs) + ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16 + ! should depend on mean radius of mode to account for gas kinetic effects + gloc = 1._kind_phys/(rhoh2o/(diff0*rhoair*qs) & + + latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair) - 1._kind_phys)) + sqrtg = sqrt(gloc) + beta = 4._kind_phys*pi*rhoh2o*gloc*gammaloc + etafactor2max = 1.e10_kind_phys/(alpha*wbar)**1.5_kind_phys ! this should make eta big if na is very small. + + do m = 1, nmode + ! skip aerosols with no dispersion, since they aren't meant to be CCN + if (dispersion_aer(m) == 0._kind_phys) then + smc(m) = 100._kind_phys + cycle + end if + ! internal mixture of aerosols + volc(m) = ma(m)/(density_aer(m)) ! only if variable size dist + if (volc(m) > 1.e-39_kind_phys .and. na(m) > 1.e-39_kind_phys) then + etafactor2(m) = 1._kind_phys/(na(m)*beta*sqrtg) !fixed or variable size dist + ! number mode radius (m) + amcubeloc(m) = (3._kind_phys*volc(m)/(4._kind_phys*pi*exp45logsig(m)*na(m))) ! only if variable size dist + smc(m) = smcrit(m) ! only for prescribed size dist + + if (hygro_aer(m) > 1.e-10_kind_phys) then ! loop only if variable size dist + smc(m) = 2._kind_phys*aten*sqrt(aten/(27._kind_phys*hygro_aer(m)*amcubeloc(m))) + else + smc(m) = 100._kind_phys + end if + else + smc(m) = 1._kind_phys + etafactor2(m) = etafactor2max ! this should make eta big if na is very small. + end if + lnsmloc(m) = log(smc(m)) ! only if variable size dist + end do + + ! single updraft + wnuc = wbar + alw = alpha*wnuc + sqrtalw = sqrt(alw) + zeta = 2._kind_phys*sqrtalw*aten/(3._kind_phys*sqrtg) + etafactor1 = 2._kind_phys*alw*sqrtalw + + do m = 1, nmode + ! skip aerosols with no dispersion, since they aren't meant to be CCN + if (dispersion_aer(m) /= 0._kind_phys) eta(m) = etafactor1*etafactor2(m) + end do + + call maxsat(zeta, eta, nmode, smc, smax) + lnsmax = log(smax) + + nact = 0._kind_phys + do m = 1, nmode + ! skip aerosols with no dispersion, since they aren't meant to be CCN + if (dispersion_aer(m) == 0._kind_phys) cycle + x = 2*(lnsmloc(m) - lnsmax)/(3*sq2*alogsig(m)) + nact = nact + 0.5_kind_phys*(1._kind_phys - erf(x))*na(m) + end do + nact = nact/rhoair ! convert from #/m3 to #/kg + + deallocate ( & + volc, & + eta, & + smc, & + etafactor2, & + amcubeloc, & + lnsmloc) + + end subroutine activate + +!=============================================================================== + + subroutine ccn_diag(ncol, pver, top_lev, maerosol, naer2, ccn_out) + use shr_spfn_mod, only: erfc => shr_spfn_erfc + + !------------------------------------------------------------------------------- + ! + ! Compute diagnostic bulk aerosol ccn concentration + ! + !------------------------------------------------------------------------------- + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev ! top level of troposphere cloud physics + real(kind_phys), intent(in) :: naer2(:, :, :) ! aerosol number concentration (1/m3) + real(kind_phys), intent(in) :: maerosol(:, :, :) ! aerosol mass conc (kg/m3) + real(kind_phys), intent(out) :: ccn_out(:, :, :) ! CCN at 6 supersaturations (#/cm3) + + ! Local variables + integer :: i, k, l, m + real(kind_phys) :: arg + + real(kind_phys) :: amcubesulfate(ncol) ! cube of dry mode radius (m) of sulfate + real(kind_phys) :: smcritsulfate(ncol) ! critical supersatuation for activation of sulfate + real(kind_phys) :: ccnfactsulfate + !------------------------------------------------------------------------------- + + ccn_out(:ncol, :, :) = 0._kind_phys + + do k = top_lev, pver + + do m = 1, naer_all + + if (m == idxsul) then + ! Lohmann treatment for sulfate has variable size distribution + do i = 1, ncol + if (naer2(i, k, m) > 0._kind_phys) then + amcubesulfate(i) = amcubefactor(m)*maerosol(i, k, m)/(naer2(i, k, m)) + smcritsulfate(i) = smcritfactor(m)/sqrt(amcubesulfate(i)) + else + smcritsulfate(i) = 1._kind_phys + end if + end do + end if + + do l = 1, psat + + if (m == idxsul) then + ! This code is modifying ccnfact for sulfate only. + do i = 1, ncol + arg = argfactor(m)*log(smcritsulfate(i)/super(l)) + if (arg < 2) then + if (arg < -2) then + ccnfactsulfate = 1.0e-6_kind_phys + else + ccnfactsulfate = 0.5e-6_kind_phys*erfc(arg) + end if + else + ccnfactsulfate = 0.0_kind_phys + end if + ccn_out(i, k, l) = ccn_out(i, k, l) + naer2(i, k, m)*ccnfactsulfate + end do + else + ! Non-sulfate species use ccnfact computed by the init routine + ccn_out(:ncol, k, l) = ccn_out(:ncol, k, l) + naer2(:ncol, k, m)*ccnfact(l, m) + end if + + end do ! supersaturation + end do ! bulk aerosol + end do ! level + + end subroutine ccn_diag + +!=============================================================================== + + subroutine maxsat(zeta, eta, nmode, smc, smax) + + ! calculates maximum supersaturation for multiple + ! competing aerosol modes. + + ! Abdul-Razzak and Ghan, A parameterization of aerosol activation. + ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844. + + real(kind_phys), intent(in) :: zeta + integer, intent(in) :: nmode ! number of modes + real(kind_phys), intent(in) :: smc(:) ! critical supersaturation for number mode radius + real(kind_phys), intent(in) :: eta(:) + + real(kind_phys), intent(out) :: smax ! maximum supersaturation + + integer :: m ! mode index + real(kind_phys) :: sum, g1, g2 + + do m = 1, nmode + if (zeta .gt. 1.e5_kind_phys*eta(m) .or. smc(m)*smc(m) .gt. 1.e5_kind_phys*eta(m)) then + ! weak forcing. essentially none activated + smax = 1.e-20_kind_phys + else + ! significant activation of this mode. calc activation all modes. + go to 1 + end if + end do + + return + +1 continue + + sum = 0 + do m = 1, nmode + if (eta(m) .gt. 1.e-20_kind_phys) then + g1 = sqrt(zeta/eta(m)) + g1 = g1*g1*g1 + g2 = smc(m)/sqrt(eta(m) + 3*zeta) + g2 = sqrt(g2) + g2 = g2*g2*g2 + sum = sum + (f1(m)*g1 + f2(m)*g2)/(smc(m)*smc(m)) + else + sum = 1.e20_kind_phys + end if + end do + + smax = 1._kind_phys/sqrt(sum) + + end subroutine maxsat + +!=============================================================================== + +end module ndrop_bam diff --git a/schemes/microp_aero/nucleate_ice.F90 b/schemes/microp_aero/nucleate_ice.F90 new file mode 100644 index 00000000..42db39a0 --- /dev/null +++ b/schemes/microp_aero/nucleate_ice.F90 @@ -0,0 +1,591 @@ +module nucleate_ice + +!------------------------------------------------------------------------------- +! Purpose: +! A parameterization of ice nucleation. +! +! *** This module is intended to be a "portable" code layer. Ideally it should +! *** not contain any use association of modules that belong to the model framework. +! +! +! Method: +! The current method is based on Liu & Penner (2005) & Liu et al. (2007) +! It related the ice nucleation with the aerosol number, temperature and the +! updraft velocity. It includes homogeneous freezing of sulfate & immersion +! freezing on mineral dust (soot disabled) in cirrus clouds, and +! Meyers et al. (1992) deposition nucleation in mixed-phase clouds +! +! The effect of preexisting ice crystals on ice nucleation in cirrus clouds is included, +! and also consider the sub-grid variability of temperature in cirrus clouds, +! following X. Shi et al. ACP (2014). +! +! Ice nucleation in mixed-phase clouds now uses classical nucleation theory (CNT), +! follows Y. Wang et al. ACP (2014), Hoose et al. (2010). +! +! Authors: +! Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010 +! Xiangjun Shi & Xiaohong Liu, 01/2014. +! +! With help from C. C. Chen and B. Eaton (2014) +!------------------------------------------------------------------------------- + +use wv_saturation, only: svp_water, svp_ice + +implicit none +private +save + +integer, parameter :: r8 = selected_real_kind(12) + +public :: nucleati_init, nucleati + +logical :: use_preexisting_ice +logical :: use_hetfrz_classnuc +logical :: use_incloud_nuc +integer :: iulog +real(r8) :: pi +real(r8) :: mincld + +real(r8), parameter :: Shet = 1.3_r8 ! het freezing threshold +real(r8), parameter :: rhoice = 0.5e3_r8 ! kg/m3, Wpice is not sensitive to rhoice +real(r8), parameter :: minweff= 0.001_r8 ! m/s +real(r8), parameter :: gamma1=1.0_r8 +real(r8), parameter :: gamma2=1.0_r8 +real(r8), parameter :: gamma3=2.0_r8 +real(r8), parameter :: gamma4=6.0_r8 + +real(r8) :: ci + +!=============================================================================== +contains +!=============================================================================== + +subroutine nucleati_init( & + use_preexisting_ice_in, use_hetfrz_classnuc_in, use_incloud_nuc_in, iulog_in, pi_in, & + mincld_in) + + logical, intent(in) :: use_preexisting_ice_in + logical, intent(in) :: use_hetfrz_classnuc_in + logical, intent(in) :: use_incloud_nuc_in + integer, intent(in) :: iulog_in + real(r8), intent(in) :: pi_in + real(r8), intent(in) :: mincld_in + + use_preexisting_ice = use_preexisting_ice_in + use_hetfrz_classnuc = use_hetfrz_classnuc_in + use_incloud_nuc = use_incloud_nuc_in + iulog = iulog_in + pi = pi_in + mincld = mincld_in + + ci = rhoice*pi/6._r8 + +end subroutine nucleati_init + +!=============================================================================== + +subroutine nucleati( & + wbar, tair, pmid, relhum, cldn, & + qc, qi, ni_in, rhoair, & + so4_num, dst_num, soot_num, subgrid, & + nuci, onihf, oniimm, onidep, onimey, & + wpice, weff, fhom, regm, & + oso4_num, odst_num, osoot_num, & + call_frm_zm_in, add_preexisting_ice_in) + + ! Input Arguments + real(r8), intent(in) :: wbar ! grid cell mean vertical velocity (m/s) + real(r8), intent(in) :: tair ! temperature (K) + real(r8), intent(in) :: pmid ! pressure at layer midpoints (pa) + real(r8), intent(in) :: relhum ! relative humidity with respective to liquid + real(r8), intent(in) :: cldn ! new value of cloud fraction (fraction) + real(r8), intent(in) :: qc ! liquid water mixing ratio (kg/kg) + real(r8), intent(in) :: qi ! grid-mean preexisting cloud ice mass mixing ratio (kg/kg) + real(r8), intent(in) :: ni_in ! grid-mean preexisting cloud ice number conc (#/kg) + real(r8), intent(in) :: rhoair ! air density (kg/m3) + real(r8), intent(in) :: so4_num ! so4 aerosol number (#/cm^3) + real(r8), intent(in) :: dst_num ! total dust aerosol number (#/cm^3) + real(r8), intent(in) :: soot_num ! soot (hydrophilic) aerosol number (#/cm^3) + real(r8), intent(in) :: subgrid ! subgrid saturation scaling factor + + ! Output Arguments + real(r8), intent(out) :: nuci ! ice number nucleated (#/kg) + real(r8), intent(out) :: onihf ! nucleated number from homogeneous freezing of so4 + real(r8), intent(out) :: oniimm ! nucleated number from immersion freezing + real(r8), intent(out) :: onidep ! nucleated number from deposition nucleation + real(r8), intent(out) :: onimey ! nucleated number from deposition nucleation (meyers: mixed phase) + real(r8), intent(out) :: wpice ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at Shom + real(r8), intent(out) :: weff ! effective Vertical velocity for ice nucleation (m/s); weff=wbar-wpice + real(r8), intent(out) :: fhom ! how much fraction of cloud can reach Shom + real(r8), intent(out) :: regm ! nucleation regime indiator + real(r8), intent(out) :: oso4_num ! so4 aerosol number (#/cm^3) + real(r8), intent(out) :: odst_num ! total dust aerosol number (#/cm^3) + real(r8), intent(out) :: osoot_num ! soot (hydrophilic) aerosol number (#/cm^3) + + ! Optional Arguments + logical, intent(in), optional :: call_frm_zm_in ! true if called from ZM convection scheme + logical, intent(in), optional :: add_preexisting_ice_in ! only false if called with pumas_v1.21+ + + ! Local workspace + real(r8) :: nihf ! nucleated number from homogeneous freezing of so4 + real(r8) :: niimm ! nucleated number from immersion freezing + real(r8) :: nidep ! nucleated number from deposition nucleation + real(r8) :: nimey ! nucleated number from deposition nucleation (meyers) + real(r8) :: n1, ni ! nucleated number + real(r8) :: tc, A, B ! work variable + real(r8) :: esl, esi, deles ! work variable + real(r8) :: wbar1, wbar2 + + ! used in SUBROUTINE Vpreice + real(r8) :: Ni_preice ! cloud ice number conc (1/m3) + real(r8) :: lami,Ri_preice ! mean cloud ice radius (m) + real(r8) :: Shom ! initial ice saturation ratio; if <1, use hom threshold Si + real(r8) :: detaT,RHimean ! temperature standard deviation, mean cloudy RHi + real(r8) :: wpicehet ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at shet + + real(r8) :: weffhet ! effective Vertical velocity for ice nucleation (m/s) weff=wbar-wpicehet + + logical :: call_frm_zm, add_preexisting_ice + !------------------------------------------------------------------------------- + + nuci = 0._r8 + + RHimean = relhum*svp_water(tair)/svp_ice(tair)*subgrid + + ! temp variables that depend on use_preexisting_ice + wbar1 = wbar + wbar2 = wbar + + ! If not using prexisting ice, the homogeneous freezing happens in the + ! entire gridbox. + fhom = 1._r8 + + if (present(call_frm_zm_in)) then + call_frm_zm = call_frm_zm_in + else + call_frm_zm = .false. + end if + + if (present(add_preexisting_ice_in)) then + add_preexisting_ice = add_preexisting_ice_in + else + add_preexisting_ice = .true. + end if + + if (use_preexisting_ice .and. (.not. call_frm_zm)) then + + Ni_preice = ni_in*rhoair ! (convert from #/kg -> #/m3) + Ni_preice = Ni_preice / max(mincld,cldn) ! in-cloud ice number density + + if (Ni_preice > 10.0_r8 .and. qi > 1.e-10_r8) then ! > 0.01/L = 10/m3 + Shom = -1.5_r8 ! if Shom<1 , Shom will be recalculated in SUBROUTINE Vpreice, according to Ren & McKenzie, 2005 + lami = (gamma4*ci*ni_in/qi)**(1._r8/3._r8) + Ri_preice = 0.5_r8/lami ! radius + Ri_preice = max(Ri_preice, 1e-8_r8) ! >0.01micron + call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shom, wpice) + call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shet, wpicehet) + else + wpice = 0.0_r8 + wpicehet = 0.0_r8 + endif + + weff = max(wbar-wpice, minweff) + wpice = min(wpice, wbar) + weffhet = max(wbar-wpicehet,minweff) + wpicehet = min(wpicehet, wbar) + + wbar1 = weff + wbar2 = weffhet + + detaT = wbar/0.23_r8 + if (use_incloud_nuc) then + call frachom(tair, 1._r8, detaT, fhom) + else + call frachom(tair, RHimean, detaT, fhom) + end if + end if + + ni = 0._r8 + tc = tair - 273.15_r8 + + ! initialize + niimm = 0._r8 + nidep = 0._r8 + nihf = 0._r8 + deles = 0._r8 + esi = 0._r8 + regm = 0._r8 + + oso4_num = 0._r8 + odst_num = 0._r8 + osoot_num = 0._r8 + + if ((so4_num >= 1.0e-10_r8 .or. (soot_num+dst_num) >= 1.0e-10_r8) .and. cldn > 0._r8) then + + if (RHimean.ge.1.2_r8) then + + if ( ((tc.le.0.0_r8).and.(tc.ge.-37.0_r8).and.(qc.lt.1.e-12_r8)).or.(tc.le.-37.0_r8)) then + + if ( (soot_num+dst_num) > 0._r8) then + A = -1.4938_r8 * log(soot_num+dst_num) + 12.884_r8 + B = -10.41_r8 * log(soot_num+dst_num) - 67.69_r8 + regm = A * log(wbar1) + B + end if + + ! heterogeneous nucleation only + if (tc .gt. regm .or. so4_num < 1.0e-10_r8) then + + if(tc.lt.-40._r8 .and. wbar1.gt.1._r8 .and. so4_num >= 1.0e-10_r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation + + call hf(tc,wbar1,relhum*subgrid,so4_num,nihf) + niimm=0._r8 + nidep=0._r8 + + ! If some homogeneous nucleation happened, assume all of the that heterogeneous + ! and coarse mode sulfate particles nucleated. + if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice + niimm = dst_num + soot_num ! assuming dst_num freeze firstly + odst_num = dst_num + osoot_num = soot_num + + oso4_num = nihf + endif + + nihf = nihf * fhom + oso4_num = oso4_num * fhom + + n1 = nihf + niimm + else + + call hetero(tc,wbar2,soot_num+dst_num,niimm,nidep) + + nihf = 0._r8 + n1 = niimm + nidep + + if ( (soot_num+dst_num) > 0._r8) then + osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num) + odst_num = dst_num * (niimm + nidep) / (soot_num + dst_num) + end if + + endif + + ! homogeneous nucleation only + else if (tc.lt.regm-5._r8 .or. (soot_num+dst_num) < 1.0e-10_r8) then + + call hf(tc,wbar1,relhum*subgrid,so4_num,nihf) + niimm=0._r8 + nidep=0._r8 + + ! If some homogeneous nucleation happened, assume all of the that + ! heterogeneous and coarse mode sulfate particles nucleated. + if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice + niimm = dst_num + soot_num ! assuming dst_num freeze firstly + odst_num = dst_num + osoot_num = soot_num + + oso4_num = nihf + endif + + nihf = nihf * fhom + oso4_num = oso4_num * fhom + + n1 = nihf + niimm + + ! transition between homogeneous and heterogeneous: interpolate in-between + else + + if (tc.lt.-40._r8 .and. wbar1.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation + + call hf(tc, wbar1, relhum*subgrid, so4_num, nihf) + niimm = 0._r8 + nidep = 0._r8 + + ! If some homogeneous nucleation happened, assume all of the + ! that heterogeneous and coarse mode sulfate particles nucleated. + if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice + niimm = dst_num + soot_num ! assuming dst_num freeze firstly + odst_num = dst_num + osoot_num = soot_num + + oso4_num = nihf + endif + + nihf = nihf * fhom + oso4_num = oso4_num * fhom + + n1 = nihf + niimm + + else + + call hf(regm-5._r8,wbar1,relhum*subgrid,so4_num,nihf) + call hetero(regm,wbar2,soot_num+dst_num,niimm,nidep) + + ! If some homogeneous nucleation happened, assume all of the + ! heterogeneous particles nucleated and add in a fraction of + ! the homogeneous freezing. + if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice + oso4_num = nihf + endif + + if ( (soot_num+dst_num) > 0._r8) then + osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num) + odst_num = dst_num * (niimm + nidep) / (soot_num + dst_num) + end if + + nihf = nihf * fhom * ((regm - tc) / 5._r8)**2 + oso4_num = oso4_num * fhom * ((regm - tc) / 5._r8)**2 + + n1 = niimm + nidep + nihf + + end if + end if + + ! Scale the rates for in-cloud number, since this is what + ! MG is expecting to find. + ni = n1 + + ! If using prexsiting ice, and allowed to add, then add it to the total. + if (use_preexisting_ice) then + if (add_preexisting_ice .and. (.not. call_frm_zm)) then + ni = ni + Ni_preice * 1e-6_r8 + end if + end if + end if + end if + end if + + ! deposition/condensation nucleation in mixed clouds (-37-64 deg) + A22_fast =-6.045_r8 !(T<=-64 deg) + B1_fast =-0.008_r8 + B21_fast =-0.042_r8 !(T>-64 deg) + B22_fast =-0.112_r8 !(T<=-64 deg) + C1_fast =0.0739_r8 + C2_fast =1.2372_r8 + + A1_slow =-0.3949_r8 + A2_slow =1.282_r8 + B1_slow =-0.0156_r8 + B2_slow =0.0111_r8 + B3_slow =0.0217_r8 + C1_slow =0.120_r8 + C2_slow =2.312_r8 + + Ni = 0.0_r8 + +!---------------------------- +!RHw parameters + A = 6.0e-4_r8*log(ww)+6.6e-3_r8 + B = 6.0e-2_r8*log(ww)+1.052_r8 + C = 1.68_r8 *log(ww)+129.35_r8 + RHw=(A*T*T+B*T+C)*0.01_r8 + + if((T.le.-37.0_r8) .and. ((RH).ge.RHw)) then + + regm = 6.07_r8*log(ww)-55.0_r8 + + if(T.ge.regm) then ! fast-growth regime + + if(T.gt.-64.0_r8) then + A2_fast=A21_fast + B2_fast=B21_fast + else + A2_fast=A22_fast + B2_fast=B22_fast + endif + + k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww)) + k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww) + + Ni = k1_fast*Na**(k2_fast) + Ni = min(Ni,Na) + + else ! slow-growth regime + + k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww)) + k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww) + + Ni = k1_slow*Na**(k2_slow) + Ni = min(Ni,Na) + + endif + + end if + +end subroutine hf + +!=============================================================================== + +SUBROUTINE Vpreice(P_in, T_in, R_in, C_in, S_in, V_out) + + ! based on Karcher et al. (2006) + ! VERTICAL VELOCITY CALCULATED FROM DEPOSITIONAL LOSS TERM + + ! SUBROUTINE arguments + REAL(r8), INTENT(in) :: P_in ! [Pa],INITIAL AIR pressure + REAL(r8), INTENT(in) :: T_in ! [K] ,INITIAL AIR temperature + REAL(r8), INTENT(in) :: R_in ! [m],INITIAL MEAN ICE CRYSTAL NUMBER RADIUS + REAL(r8), INTENT(in) :: C_in ! [m-3],INITIAL TOTAL ICE CRYSTAL NUMBER DENSITY, [1/cm3] + REAL(r8), INTENT(in) :: S_in ! [-],INITIAL ICE SATURATION RATIO;; if <1, use hom threshold Si + REAL(r8), INTENT(out) :: V_out ! [m/s], VERTICAL VELOCITY REDUCTION (caused by preexisting ice) + + ! SUBROUTINE parameters + REAL(r8), PARAMETER :: ALPHAc = 0.5_r8 ! density of ice (g/cm3), !!!V is not related to ALPHAc + REAL(r8), PARAMETER :: FA1c = 0.601272523_r8 + REAL(r8), PARAMETER :: FA2c = 0.000342181855_r8 + REAL(r8), PARAMETER :: FA3c = 1.49236645E-12_r8 + REAL(r8), PARAMETER :: WVP1c = 3.6E+10_r8 + REAL(r8), PARAMETER :: WVP2c = 6145.0_r8 + REAL(r8), PARAMETER :: FVTHc = 11713803.0_r8 + REAL(r8), PARAMETER :: THOUBKc = 7.24637701E+18_r8 + REAL(r8), PARAMETER :: SVOLc = 3.23E-23_r8 ! SVOL=XMW/RHOICE + REAL(r8), PARAMETER :: FDc = 249.239822_r8 + REAL(r8), PARAMETER :: FPIVOLc = 3.89051704E+23_r8 + REAL(r8) :: T,P,S,R,C + REAL(r8) :: A1,A2,A3,B1,B2 + REAL(r8) :: T_1,PICE,FLUX,ALP4,CISAT,DLOSS,VICE + + T = T_in ! K , K + P = P_in*1e-2_r8 ! Pa , hpa + + IF (S_in.LT.1.0_r8) THEN + S = 2.349_r8 - (T/259.0_r8) ! homogeneous freezing threshold, according to Ren & McKenzie, 2005 + ELSE + S = S_in ! INPUT ICE SATURATION RATIO, -, >1 + ENDIF + + R = R_in*1e2_r8 ! m => cm + C = C_in*1e-6_r8 ! m-3 => cm-3 + T_1 = 1.0_r8/ T + PICE = WVP1c * EXP(-(WVP2c*T_1)) + ALP4 = 0.25_r8 * ALPHAc + FLUX = ALP4 * SQRT(FVTHc*T) + CISAT = THOUBKc * PICE * T_1 + A1 = ( FA1c * T_1 - FA2c ) * T_1 + A2 = 1.0_r8/ CISAT + A3 = FA3c * T_1 / P + B1 = FLUX * SVOLc * CISAT * ( S-1.0_r8 ) + B2 = FLUX * FDc * P * T_1**1.94_r8 + DLOSS = FPIVOLc * C * B1 * R**2 / ( 1.0_r8+ B2 * R ) + VICE = ( A2 + A3 * S ) * DLOSS / ( A1 * S ) ! 2006,(19) + V_out = VICE*1e-2_r8 ! cm/s => m/s + +END SUBROUTINE Vpreice + +subroutine frachom(Tmean,RHimean,detaT,fhom) + ! How much fraction of cirrus might reach Shom + ! base on "A cirrus cloud scheme for general circulation models", + ! B. Karcher and U. Burkhardt 2008 + + real(r8), intent(in) :: Tmean, RHimean, detaT + real(r8), intent(out) :: fhom + + real(r8), parameter :: seta = 6132.9_r8 ! K + integer, parameter :: Nbin=200 ! (Tmean - 3*detaT, Tmean + 3*detaT) + + real(r8) :: PDF_T(Nbin) ! temperature PDF; ! PDF_T=0 outside (Tmean-3*detaT, Tmean+3*detaT) + real(r8) :: Sbin(Nbin) ! the fluctuations of Si that are driven by the T variations + real(r8) :: Sihom, deta + integer :: i + + Sihom = 2.349_r8-Tmean/259.0_r8 ! homogeneous freezing threshold, according to Ren & McKenzie, 2005 + fhom = 0.0_r8 + + do i = Nbin, 1, -1 + + deta = (i - 0.5_r8 - Nbin/2)*6.0_r8/Nbin ! PDF_T=0 outside (Tmean-3*detaT, Tmean+3*detaT) + Sbin(i) = RHimean*exp(deta*detaT*seta/Tmean**2.0_r8) + PDF_T(i) = exp(-deta**2.0_r8/2.0_r8)*6.0_r8/(sqrt(2.0_r8*Pi)*Nbin) + + + if (Sbin(i).ge.Sihom) then + fhom = fhom + PDF_T(i) + else + exit + end if + end do + + fhom = min(1.0_r8, fhom/0.997_r8) ! accounting for the finite limits (-3 , 3) +end subroutine frachom + +end module nucleate_ice diff --git a/schemes/microp_aero/scale_subgrid_vertical_velocity.F90 b/schemes/microp_aero/scale_subgrid_vertical_velocity.F90 new file mode 100644 index 00000000..42f82227 --- /dev/null +++ b/schemes/microp_aero/scale_subgrid_vertical_velocity.F90 @@ -0,0 +1,68 @@ +! Apply min/max/scale to subgrid vertical velocity and derive wsubi. +module scale_subgrid_vertical_velocity + + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: scale_subgrid_vertical_velocity_run + +contains + +!> \section arg_table_scale_subgrid_vertical_velocity_run Argument Table +!! \htmlinclude scale_subgrid_vertical_velocity_run.html + subroutine scale_subgrid_vertical_velocity_run( & + ncol, pver, top_lev, & + wsub_min, wsubi_min, & + wsub_scale, wsubi_scale, & + use_preexisting_ice, & + wsub, wsubi, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev ! top vertical level for cloud physics + + real(kind_phys), intent(in) :: wsub_min ! minimum subgrid vertical velocity for droplet nucleation [m s-1] + real(kind_phys), intent(in) :: wsubi_min ! minimum subgrid vertical velocity for ice nucleation [m s-1] + real(kind_phys), intent(in) :: wsub_scale ! scaling factor for subgrid vertical velocity (liquid) + real(kind_phys), intent(in) :: wsubi_scale ! scaling factor for subgrid vertical velocity (ice) + + logical, intent(in) :: use_preexisting_ice ! use preexisting ice? [flag] + + ! Input/Output arguments + real(kind_phys), intent(inout) :: wsub(:, :) ! subgrid vertical velocity [m s-1] (ncol, pver) + + ! Output arguments + real(kind_phys), intent(out) :: wsubi(:, :) ! subgrid vertical velocity for ice nucleation [m s-1] (ncol, pver) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + + errmsg = '' + errflg = 0 + + ! Set minimum values above top_lev + wsub (:, :top_lev-1) = wsub_min + wsubi(:, :top_lev-1) = wsubi_min + + do k = top_lev, pver + do i = 1, ncol + ! Derive wsubi from raw wsub BEFORE scaling wsub + wsubi(i,k) = max(wsubi_min, wsub(i,k)) * wsubi_scale + if (.not. use_preexisting_ice) then + wsubi(i,k) = min(wsubi(i,k), 0.2_kind_phys) + end if + + wsub(i,k) = max(wsub_min, wsub(i,k)) * wsub_scale + end do + end do + + end subroutine scale_subgrid_vertical_velocity_run + +end module scale_subgrid_vertical_velocity diff --git a/schemes/microp_aero/scale_subgrid_vertical_velocity.meta b/schemes/microp_aero/scale_subgrid_vertical_velocity.meta new file mode 100644 index 00000000..0f37142a --- /dev/null +++ b/schemes/microp_aero/scale_subgrid_vertical_velocity.meta @@ -0,0 +1,79 @@ +[ccpp-table-properties] + name = scale_subgrid_vertical_velocity + type = scheme + +[ccpp-arg-table] + name = scale_subgrid_vertical_velocity_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ wsub_min ] + standard_name = minimum_subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ wsubi_min ] + standard_name = minimum_subgrid_scale_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ wsub_scale ] + standard_name = tunable_parameter_for_subgrid_vertical_velocity_scaling_for_droplet_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ wsubi_scale ] + standard_name = tunable_parameter_for_subgrid_vertical_velocity_scaling_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ use_preexisting_ice ] + standard_name = do_preexisting_ice_treatment_in_ice_nucleation + units = flag + type = logical + dimensions = () + intent = in +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = inout +[ wsubi ] + standard_name = subgrid_scale_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/microp_aero/scale_subgrid_vertical_velocity_namelist.xml b/schemes/microp_aero/scale_subgrid_vertical_velocity_namelist.xml new file mode 100644 index 00000000..9e582e52 --- /dev/null +++ b/schemes/microp_aero/scale_subgrid_vertical_velocity_namelist.xml @@ -0,0 +1,73 @@ + + + + + + + + real + microphys + wsub_scale_nl + minimum_subgrid_scale_vertical_velocity_for_droplet_nucleation + m s-1 + + Minimum subgrid vertical velocity (before scale factor) for liquid droplet + nucleation during aerosol activation. + Default: 0.2 (CAM6), 0.1 (CAM7) + + + 0.2D0 + 0.1D0 + + + + + real + microphys + wsub_scale_nl + minimum_subgrid_scale_vertical_velocity_for_ice_nucleation + m s-1 + + Minimum subgrid vertical velocity for ice nucleation during aerosol activation. + Default: 0.001 + + + 0.001D0 + + + + + real + microphys + wsub_scale_nl + tunable_parameter_for_subgrid_vertical_velocity_scaling_for_droplet_nucleation + 1 + + Unitless scaling factor for the liquid droplet subgrid scale vertical velocity + during aerosol activation. + Default: 1.0 (CAM6), 0.9 (CAM7) + + + 1.D0 + 0.9D0 + + + + + real + microphys + wsub_scale_nl + tunable_parameter_for_subgrid_vertical_velocity_scaling_for_ice_nucleation + 1 + + Unitless scaling factor for the ice droplet subgrid scale vertical velocity + during aerosol activation. + Default: 1.0 (CAM6), 1.5 (CAM7) + + + 1.D0 + 1.5D0 + + + + diff --git a/schemes/rrtmgp/rrtmgp_lw_aerosols.F90 b/schemes/rrtmgp/rrtmgp_lw_aerosols.F90 index 4078204a..7ab54daf 100644 --- a/schemes/rrtmgp/rrtmgp_lw_aerosols.F90 +++ b/schemes/rrtmgp/rrtmgp_lw_aerosols.F90 @@ -12,18 +12,27 @@ module rrtmgp_lw_aerosols !> \section arg_table_rrtmgp_lw_aerosols_run Argument Table !! \htmlinclude rrtmgp_lw_aerosols_run.html !! - subroutine rrtmgp_lw_aerosols_run(dolwrad, aer_lw, errmsg, errflg) + subroutine rrtmgp_lw_aerosols_run( & + dolwrad, ncol, ktopcam, ktoprad, & + aer_lw_abs, & + aer_lw, errmsg, errflg) + use ccpp_optical_props, only: ty_optical_props_1scl_ccpp use ccpp_kinds, only: kind_phys ! Inputs - logical, intent(in) :: dolwrad !< Flag to perform longwave calculation + logical, intent(in) :: dolwrad !< do longwave calculation + integer, intent(in) :: ncol !< number of columns + integer, intent(in) :: ktopcam !< top CAM level for RRTMGP + integer, intent(in) :: ktoprad !< top RRTMGP level + + real(kind_phys), intent(in) :: aer_lw_abs(:, :, :) !< LW absorption OD (ncol,pver,nlwbands) ! Outputs - class(ty_optical_props_1scl_ccpp), intent(inout) :: aer_lw !< Aerosol optical properties object + class(ty_optical_props_1scl_ccpp), intent(inout) :: aer_lw !< Aerosol optical properties object - character(len=*), intent(out) :: errmsg !< CCPP error message - integer, intent(out) :: errflg !< CCPP error flag + character(len=*), intent(out) :: errmsg !< CCPP error message + integer, intent(out) :: errflg !< CCPP error flag ! Initialize CCPP error handling variables errmsg = '' @@ -31,9 +40,12 @@ subroutine rrtmgp_lw_aerosols_run(dolwrad, aer_lw, errmsg, errflg) if (.not. dolwrad) return - ! REMOVECAM: scheme is just a stub for now - ! Set the optical properties to 0 for snapshot testing + ! Initialize to zero (provides zero optical depth in any extra layer) aer_lw%optical_props%tau = 0.0_kind_phys + ! Map CAM vertical levels to RRTMGP vertical levels + ! No band reorder needed for LW + aer_lw%optical_props%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :) + end subroutine rrtmgp_lw_aerosols_run end module rrtmgp_lw_aerosols diff --git a/schemes/rrtmgp/rrtmgp_lw_aerosols.meta b/schemes/rrtmgp/rrtmgp_lw_aerosols.meta index f5598619..f874cf05 100644 --- a/schemes/rrtmgp/rrtmgp_lw_aerosols.meta +++ b/schemes/rrtmgp/rrtmgp_lw_aerosols.meta @@ -11,6 +11,30 @@ type = logical dimensions = () intent = in +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ ktopcam ] + standard_name = vertical_index_at_top_level_where_RRTMGP_is_active + units = index + type = integer + dimensions = () + intent = in +[ ktoprad ] + standard_name = RRTMGP_vertical_index_at_top_level_where_RRTMGP_is_active + units = index + type = integer + dimensions = () + intent = in +[ aer_lw_abs ] + standard_name = longwave_aerosol_absorption_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_longwave_radiation) + intent = in [ aer_lw ] standard_name = longwave_aerosol_optical_properties_object_for_RRTMGP units = none diff --git a/schemes/rrtmgp/rrtmgp_sw_aerosols.F90 b/schemes/rrtmgp/rrtmgp_sw_aerosols.F90 index 5f72395c..9f1f49d4 100644 --- a/schemes/rrtmgp/rrtmgp_sw_aerosols.F90 +++ b/schemes/rrtmgp/rrtmgp_sw_aerosols.F90 @@ -3,28 +3,57 @@ !> This module sets the RRTMGP aerosol shortwave optical properties module rrtmgp_sw_aerosols + + use ccpp_kinds, only: kind_phys + implicit none private public rrtmgp_sw_aerosols_run + + ! Mapping from RRTMG shortwave bands to RRTMGP + integer, parameter, dimension(14) :: rrtmg_to_rrtmgp_swbands = & + [ 14, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ] + + real(kind_phys), parameter :: tiny = 1.0e-80_kind_phys + contains !> \section arg_table_rrtmgp_sw_aerosols_run Argument Table !! \htmlinclude rrtmgp_sw_aerosols_run.html !! - subroutine rrtmgp_sw_aerosols_run(doswrad, nday, aer_sw, errmsg, errflg) + subroutine rrtmgp_sw_aerosols_run( & + doswrad, ncol, pver, nswbands, nday, idxday, ktopcam, ktoprad, & + aer_tau, aer_tau_w, aer_tau_w_g, & + aer_sw, errmsg, errflg) + use ccpp_optical_props, only: ty_optical_props_2str_ccpp - use ccpp_kinds, only: kind_phys ! Inputs - logical, intent(in) :: doswrad !< Flag to perform shortwave calculation - integer, intent(in) :: nday !< daytime points dimension + logical, intent(in) :: doswrad !< do shortwave calculation + integer, intent(in) :: ncol !< number of columns + integer, intent(in) :: pver !< number of vertical layers + integer, intent(in) :: nswbands !< number of SW bands + integer, intent(in) :: nday !< daytime columns count + integer, intent(in) :: idxday(:) !< daytime column indices + integer, intent(in) :: ktopcam !< top CAM level for RRTMGP + integer, intent(in) :: ktoprad !< top RRTMGP level + + real(kind_phys), intent(in) :: aer_tau(:, :, :) !< SW extinction OD (ncol,pver,nswbands) + real(kind_phys), intent(in) :: aer_tau_w(:, :, :) !< SW ssa*tau (ncol,pver,nswbands) + real(kind_phys), intent(in) :: aer_tau_w_g(:, :, :) !< SW asy*ssa*tau (ncol,pver,nswbands) ! Outputs - class(ty_optical_props_2str_ccpp), intent(inout) :: aer_sw !< Aerosol optical properties object + class(ty_optical_props_2str_ccpp), intent(inout) :: aer_sw !< Aerosol optical properties object - character(len=*), intent(out) :: errmsg !< CCPP error message - integer, intent(out) :: errflg !< CCPP error flag + character(len=*), intent(out) :: errmsg !< CCPP error message + integer, intent(out) :: errflg !< CCPP error flag + + ! Local variables + integer :: i + real(kind_phys) :: tau_reordered(ncol, pver, nswbands) + real(kind_phys) :: tau_w_reordered(ncol, pver, nswbands) + real(kind_phys) :: tau_w_g_reordered(ncol, pver, nswbands) ! Initialize CCPP error handling variables errmsg = '' @@ -32,11 +61,44 @@ subroutine rrtmgp_sw_aerosols_run(doswrad, nday, aer_sw, errmsg, errflg) if (.not. doswrad .or. nday == 0) return - ! REMOVECAM: scheme is just a stub for now - ! Set the optical properties to constants for snapshot testing + ! Reorder from RRTMG band order to RRTMGP band order + tau_reordered(:ncol, :, :) = aer_tau(:ncol, :, rrtmg_to_rrtmgp_swbands) + tau_w_reordered(:ncol, :, :) = aer_tau_w(:ncol, :, rrtmg_to_rrtmgp_swbands) + tau_w_g_reordered(:ncol, :, :) = aer_tau_w_g(:ncol, :, rrtmg_to_rrtmgp_swbands) + + ! Initialize defaults: tau=0, ssa=1, g=0 aer_sw%optical_props%tau = 0.0_kind_phys - aer_sw%optical_props%g = 0.0_kind_phys aer_sw%optical_props%ssa = 1.0_kind_phys + aer_sw%optical_props%g = 0.0_kind_phys + + ! Subset to daytime columns and map CAM vertical to RRTMGP vertical + do i = 1, nday + ! Set aerosol optical depth, clip to zero + aer_sw%optical_props%tau(i, ktoprad:, :) = & + max(tau_reordered(idxday(i), ktopcam:, :), 0._kind_phys) + + ! Set single scattering albedo: ssa = tau_w / tau + where (tau_reordered(idxday(i), ktopcam:, :) > 0._kind_phys) + aer_sw%optical_props%ssa(i, ktoprad:, :) = & + tau_w_reordered(idxday(i), ktopcam:, :) & + / tau_reordered(idxday(i), ktopcam:, :) + elsewhere + aer_sw%optical_props%ssa(i, ktoprad:, :) = 1._kind_phys + end where + + ! Set asymmetry parameter: g = tau_w_g / tau_w + where (tau_w_reordered(idxday(i), ktopcam:, :) > tiny) + aer_sw%optical_props%g(i, ktoprad:, :) = & + tau_w_g_reordered(idxday(i), ktopcam:, :) & + / tau_w_reordered(idxday(i), ktopcam:, :) + elsewhere + aer_sw%optical_props%g(i, ktoprad:, :) = 0._kind_phys + end where + end do + + ! Impose limits on the components + aer_sw%optical_props%ssa = min(max(aer_sw%optical_props%ssa, 0._kind_phys), 1._kind_phys) + aer_sw%optical_props%g = min(max(aer_sw%optical_props%g, -1._kind_phys), 1._kind_phys) end subroutine rrtmgp_sw_aerosols_run end module rrtmgp_sw_aerosols diff --git a/schemes/rrtmgp/rrtmgp_sw_aerosols.meta b/schemes/rrtmgp/rrtmgp_sw_aerosols.meta index f9736f22..2eff769f 100644 --- a/schemes/rrtmgp/rrtmgp_sw_aerosols.meta +++ b/schemes/rrtmgp/rrtmgp_sw_aerosols.meta @@ -11,12 +11,66 @@ type = logical dimensions = () intent = in +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ nswbands ] + standard_name = number_of_bands_for_shortwave_radiation + units = count + type = integer + dimensions = () + intent = in [ nday ] standard_name = daytime_columns_dimension units = count type = integer dimensions = () intent = in +[ idxday ] + standard_name = daytime_columns + units = index + type = integer + dimensions = (horizontal_loop_extent) + intent = in +[ ktopcam ] + standard_name = vertical_index_at_top_level_where_RRTMGP_is_active + units = index + type = integer + dimensions = () + intent = in +[ ktoprad ] + standard_name = RRTMGP_vertical_index_at_top_level_where_RRTMGP_is_active + units = index + type = integer + dimensions = () + intent = in +[ aer_tau ] + standard_name = shortwave_aerosol_extinction_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = in +[ aer_tau_w ] + standard_name = shortwave_aerosol_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = in +[ aer_tau_w_g ] + standard_name = shortwave_aerosol_asymmetry_parameter_times_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = in [ aer_sw ] standard_name = shortwave_aerosol_optical_properties_object_for_RRTMGP units = none diff --git a/schemes/sima_diagnostics/aerosol_optics_diagnostics.F90 b/schemes/sima_diagnostics/aerosol_optics_diagnostics.F90 new file mode 100644 index 00000000..79ffc797 --- /dev/null +++ b/schemes/sima_diagnostics/aerosol_optics_diagnostics.F90 @@ -0,0 +1,458 @@ +! Aerosol optics diagnostics scheme for CAM-SIMA. +! Outputs aggregate SW/LW optical property diagnostics and per-species +! AODs/burdens from the aerosol_optics CCPP scheme. +! +! Day-only fields use fillvalue for night columns; +! day+night (dn) fields carry data for all columns. +! +! Author: Haipeng Lin, NSF-NCAR/CGD/AMP, March 2026 +module aerosol_optics_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: aerosol_optics_diagnostics_init + public :: aerosol_optics_diagnostics_run + + ! ODV_ aerosol diagnostics for bulk aerosol + integer :: num_odv = 0 + character(len=64), allocatable :: odv_names(:) + +contains + + !> \section arg_table_aerosol_optics_diagnostics_init Argument Table + !! \htmlinclude aerosol_optics_diagnostics_init.html + subroutine aerosol_optics_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + use cam_history_support, only: horiz_only + use radiative_aerosol_definitions, only: bulk_aerosol_list + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: i + + errmsg = '' + errflg = 0 + + !----------------------------------------------------------------- + ! SW column diagnostics (horiz_only, avg) + !----------------------------------------------------------------- + call history_add_field('AODVIS', 'Aerosol optical depth at 550 nm', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODVISdn', 'Aerosol optical depth at 550 nm, day and night', horiz_only, 'avg', '1') + call history_add_field('AODUV', 'Aerosol optical depth in UV band', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODUVdn', 'Aerosol optical depth in UV band, day and night', horiz_only, 'avg', '1') + call history_add_field('AODNIR', 'Aerosol optical depth in near-IR band', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODNIRdn', 'Aerosol optical depth in near-IR band, day and night', horiz_only, 'avg', '1') + call history_add_field('AODABS', 'Aerosol absorption optical depth at 550 nm', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODABSdn', 'Aerosol absorption optical depth at 550 nm, day and night', horiz_only, 'avg', '1') + call history_add_field('AODTOT', 'Total aerosol optical depth (all SW bands)', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODTOTdn', 'Total aerosol optical depth (all SW bands), day and night', horiz_only, 'avg', '1') + call history_add_field('SSAVIS', 'Aerosol single scatter albedo at 550 nm', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('SSAVISdn', 'Aerosol single scatter albedo at 550 nm, day and night', horiz_only, 'avg', '1') + + ! Stratospheric column AODs + call history_add_field('AODVISst', 'Stratospheric aerosol optical depth at 550 nm', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODVISstdn', 'Stratospheric aerosol optical depth at 550 nm, day and night', horiz_only, 'avg', '1') + call history_add_field('AODUVst', 'Stratospheric aerosol optical depth in UV band', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODUVstdn', 'Stratospheric aerosol optical depth in UV band, day and night', horiz_only, 'avg', '1') + call history_add_field('AODNIRst', 'Stratospheric aerosol optical depth in near-IR band', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODNIRstdn', 'Stratospheric aerosol optical depth in near-IR band, day and night', horiz_only, 'avg', '1') + + !----------------------------------------------------------------- + ! SW per-level diagnostics (lev, avg) + !----------------------------------------------------------------- + call history_add_field('EXTINCT', 'Aerosol extinction at 550 nm', 'lev', 'avg', 'm-1', flag_xyfill=.true.) + call history_add_field('EXTINCTdn', 'Aerosol extinction at 550 nm, day and night', 'lev', 'avg', 'm-1') + call history_add_field('EXTINCTUV', 'Aerosol extinction in UV band', 'lev', 'avg', 'm-1', flag_xyfill=.true.) + call history_add_field('EXTINCTUVdn', 'Aerosol extinction in UV band, day and night', 'lev', 'avg', 'm-1') + call history_add_field('EXTINCTNIR', 'Aerosol extinction in near-IR band', 'lev', 'avg', 'm-1', flag_xyfill=.true.) + call history_add_field('EXTINCTNIRdn', 'Aerosol extinction in near-IR band, day and night', 'lev', 'avg', 'm-1') + call history_add_field('ABSORB', 'Aerosol absorption at 550 nm', 'lev', 'avg', 'm-1', flag_xyfill=.true.) + call history_add_field('ABSORBdn', 'Aerosol absorption at 550 nm, day and night', 'lev', 'avg', 'm-1') + + !----------------------------------------------------------------- + ! Per-species column diagnostics (horiz_only, avg) + !----------------------------------------------------------------- + call history_add_field('AODDUST', 'Dust aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODDUSTdn', 'Dust aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODSO4', 'Sulfate aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODSO4dn', 'Sulfate aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODBC', 'Black carbon aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODBCdn', 'Black carbon aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODPOM', 'POM aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODPOMdn', 'POM aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODSOA', 'SOA aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODSOAdn', 'SOA aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODSS', 'Seasalt aerosol optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODSSdn', 'Seasalt aerosol optical depth, day and night', horiz_only, 'avg', '1') + call history_add_field('AODABSBC', 'Black carbon absorption optical depth', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODABSBCdn', 'Black carbon absorption optical depth, day and night', horiz_only, 'avg', '1') + + call history_add_field('BURDENDUST', 'Dust column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENDUSTdn', 'Dust column burden, day and night', horiz_only, 'avg', 'kg m-2') + call history_add_field('BURDENSO4', 'Sulfate column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENSO4dn', 'Sulfate column burden, day and night', horiz_only, 'avg', 'kg m-2') + call history_add_field('BURDENBC', 'Black carbon column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENBCdn', 'Black carbon column burden, day and night', horiz_only, 'avg', 'kg m-2') + call history_add_field('BURDENPOM', 'POM column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENPOMdn', 'POM column burden, day and night', horiz_only, 'avg', 'kg m-2') + call history_add_field('BURDENSOA', 'SOA column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENSOAdn', 'SOA column burden, day and night', horiz_only, 'avg', 'kg m-2') + call history_add_field('BURDENSEASALT', 'Seasalt column burden', horiz_only, 'avg', 'kg m-2', flag_xyfill=.true.) + call history_add_field('BURDENSEASALTdn','Seasalt column burden, day and night', horiz_only, 'avg', 'kg m-2') + + !----------------------------------------------------------------- + ! LW per-level diagnostics (lev, avg) + !----------------------------------------------------------------- + call history_add_field('TOTABSLW', 'Total LW aerosol absorption optical depth', 'lev', 'avg', '1') + call history_add_field('AODABSLW', 'LW aerosol absorption optical depth at diagnostic band', 'lev', 'avg', '1') + + !----------------------------------------------------------------- + ! Additional CAM-standard column diagnostics (aliases) + !----------------------------------------------------------------- + call history_add_field('AEROD_v', 'Total aerosol optical depth in visible band', horiz_only, 'avg', '1', flag_xyfill=.true.) + call history_add_field('AODvstrt', 'Stratospheric aerosol optical depth in visible band', horiz_only, 'avg', '1', flag_xyfill=.true.) + + !----------------------------------------------------------------- + ! Per-aerosol visible OD diagnostics (ODV_) for bulk aerosols + !----------------------------------------------------------------- + num_odv = bulk_aerosol_list(0)%numaerosols + if (num_odv > 0) then + allocate(odv_names(num_odv)) + do i = 1, num_odv + odv_names(i) = 'ODV_' // trim(bulk_aerosol_list(0)%aer(i)%camname) + call history_add_field(trim(odv_names(i)), & + trim(bulk_aerosol_list(0)%aer(i)%camname)//' optical depth in visible band', & + horiz_only, 'avg', '1', flag_xyfill=.true.) + end do + end if + + end subroutine aerosol_optics_diagnostics_init + + !> \section arg_table_aerosol_optics_diagnostics_run Argument Table + !! \htmlinclude aerosol_optics_diagnostics_run.html + subroutine aerosol_optics_diagnostics_run( & + ncol, pver, nswbands, nlwbands, & + N_DIAG, active_calls, & + constituents, & + aer_tau, aer_tau_w, aer_lw_abs, & + idx_sw_diag, idx_uv_diag, idx_nir_diag, idx_lw_diag, & + pdeldry, pmid, t, rga, rair, & + nnite, idxnite, troplev, & + dustaod, sulfaod, bcaod, pomaod, soaaod, ssltaod, aodabsbc, & + burdendust, burdenso4, burdenbc, burdenpom, burdensoa, burdenseasalt, & + ssavis, aodvis, & + num_bulk_aer, odv_col_aod, & + errmsg, errflg) + + use cam_history, only: history_out_field + use cam_history_support, only: fillvalue + + ! host-dependent diagnostic output subroutine for radiatively active aerosol: + use aerosol_mmr_host, only: rad_aer_diag_out + + !----------------------------------------------------------------- + ! Input arguments + !----------------------------------------------------------------- + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: nswbands + integer, intent(in) :: nlwbands + integer, intent(in) :: N_DIAG ! max number of diagnostic lists + logical, intent(in) :: active_calls(0:) ! flag for active diagnostic list calls + ! n.b.: the lower bound of the active_calls array has to be set here. + ! in radiative_aerosol_definitions, active_calls is 0-indexed because + ! 0 is the climate list index. + ! The Fortran 90 standard (and higher) states assumed-shape arrays + ! lower bound default to 1 unless specified (ISO/IEC 1539 5.1.2.4.2) + + real(kind_phys), intent(in) :: constituents(:, :, :) + + real(kind_phys), intent(in) :: aer_tau(:,:,:) ! SW extinction OD (ncol,pver,nswbands) + real(kind_phys), intent(in) :: aer_tau_w(:,:,:) ! SW ssa*tau (ncol,pver,nswbands) + real(kind_phys), intent(in) :: aer_lw_abs(:,:,:) ! LW absorption OD (ncol,pver,nlwbands) + integer, intent(in) :: idx_sw_diag ! Index of SW diagnostic (vis) band + integer, intent(in) :: idx_uv_diag ! Index of UV band + integer, intent(in) :: idx_nir_diag ! Index of near-IR band + integer, intent(in) :: idx_lw_diag ! Index of LW diagnostic band + real(kind_phys), intent(in) :: pdeldry(:,:) ! Dry air pressure thickness (ncol,pver) [Pa] + real(kind_phys), intent(in) :: pmid(:,:) ! Air pressure (ncol,pver) [Pa] + real(kind_phys), intent(in) :: t(:,:) ! Air temperature (ncol,pver) [K] + real(kind_phys), intent(in) :: rga ! 1/g [s2 m-1] + real(kind_phys), intent(in) :: rair ! Gas constant of dry air [J kg-1 K-1] + integer, intent(in) :: nnite ! Number of night columns + integer, intent(in) :: idxnite(:) ! Indices of night columns + integer, intent(in) :: troplev(:) ! Tropopause level index (ncol) + + ! Per-species inputs from aerosol_optics + real(kind_phys), intent(in) :: dustaod(:) + real(kind_phys), intent(in) :: sulfaod(:) + real(kind_phys), intent(in) :: bcaod(:) + real(kind_phys), intent(in) :: pomaod(:) + real(kind_phys), intent(in) :: soaaod(:) + real(kind_phys), intent(in) :: ssltaod(:) + real(kind_phys), intent(in) :: aodabsbc(:) + real(kind_phys), intent(in) :: burdendust(:) + real(kind_phys), intent(in) :: burdenso4(:) + real(kind_phys), intent(in) :: burdenbc(:) + real(kind_phys), intent(in) :: burdenpom(:) + real(kind_phys), intent(in) :: burdensoa(:) + real(kind_phys), intent(in) :: burdenseasalt(:) + real(kind_phys), intent(in) :: ssavis(:) + real(kind_phys), intent(in) :: aodvis(:) + + ! Per-constituent visible OD inputs from aerosol_optics + integer, intent(in) :: num_bulk_aer + real(kind_phys), intent(in) :: odv_col_aod(:,:) ! (ncol, num_bulk_aer) + + ! CCPP error handling + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + integer :: icol, i, iaer + integer :: ilist + + ! Derived column quantities + real(kind_phys) :: aodvis_l(ncol) + real(kind_phys) :: aoduv(ncol) + real(kind_phys) :: aodnir(ncol) + real(kind_phys) :: aodabs(ncol) + real(kind_phys) :: aodtot(ncol) + real(kind_phys) :: ssavis_l(ncol) + + ! Stratospheric column AODs + real(kind_phys) :: aodvisst(ncol) + real(kind_phys) :: aoduvst(ncol) + real(kind_phys) :: aodnirst(ncol) + + ! Per-level quantities + real(kind_phys) :: mass(ncol, pver) + real(kind_phys) :: air_density(ncol, pver) + real(kind_phys) :: extinct(ncol, pver) + real(kind_phys) :: extinctuv(ncol, pver) + real(kind_phys) :: extinctnir(ncol, pver) + real(kind_phys) :: absorb(ncol, pver) + + ! LW diagnostics + real(kind_phys) :: totabslw(ncol, pver) + real(kind_phys) :: aodabslw(ncol, pver) + + ! Local copy for ODV_ output with fillvalue + real(kind_phys) :: odv_tmp(ncol) + + ! Local copies of per-species fields (for fillvalue modification) + real(kind_phys) :: dustaod_l(ncol), sulfaod_l(ncol), bcaod_l(ncol) + real(kind_phys) :: pomaod_l(ncol), soaaod_l(ncol), ssltaod_l(ncol) + real(kind_phys) :: aodabsbc_l(ncol) + real(kind_phys) :: burdendust_l(ncol), burdenso4_l(ncol), burdenbc_l(ncol) + real(kind_phys) :: burdenpom_l(ncol), burdensoa_l(ncol), burdenseasalt_l(ncol) + + errmsg = '' + errflg = 0 + + !----------------------------------------------------------------- + ! Derived quantities + !----------------------------------------------------------------- + mass(:ncol, :) = pdeldry(:ncol, :) * rga + air_density(:ncol, :) = pmid(:ncol, :) / (rair * t(:ncol, :)) + + !----------------------------------------------------------------- + ! Column AODs (sum over levels) + !----------------------------------------------------------------- + ! aodvis and ssavis are already computed in aerosol_optics_run + aodvis_l(:ncol) = aodvis(:ncol) + ssavis_l(:ncol) = ssavis(:ncol) + + do icol = 1, ncol + aoduv(icol) = sum(aer_tau(icol, :, idx_uv_diag)) + aodnir(icol) = sum(aer_tau(icol, :, idx_nir_diag)) + aodabs(icol) = sum(aer_tau(icol, :, idx_sw_diag) - aer_tau_w(icol, :, idx_sw_diag)) + aodtot(icol) = sum(aer_tau(icol, :, :)) + end do + + !----------------------------------------------------------------- + ! Stratospheric column AODs (levels above tropopause) + !----------------------------------------------------------------- + do icol = 1, ncol + aodvisst(icol) = sum(aer_tau(icol, 1:troplev(icol), idx_sw_diag)) + aoduvst(icol) = sum(aer_tau(icol, 1:troplev(icol), idx_uv_diag)) + aodnirst(icol) = sum(aer_tau(icol, 1:troplev(icol), idx_nir_diag)) + end do + + !----------------------------------------------------------------- + ! Per-level extinction [m-1] and absorption [m-1] + !----------------------------------------------------------------- + extinct(:ncol, :) = aer_tau(:ncol, :, idx_sw_diag) * air_density(:ncol, :) / mass(:ncol, :) + extinctuv(:ncol, :) = aer_tau(:ncol, :, idx_uv_diag) * air_density(:ncol, :) / mass(:ncol, :) + extinctnir(:ncol, :) = aer_tau(:ncol, :, idx_nir_diag) * air_density(:ncol, :) / mass(:ncol, :) + absorb(:ncol, :) = (aer_tau(:ncol, :, idx_sw_diag) - aer_tau_w(:ncol, :, idx_sw_diag)) & + * air_density(:ncol, :) / mass(:ncol, :) + + !----------------------------------------------------------------- + ! LW diagnostics + !----------------------------------------------------------------- + totabslw(:ncol, :) = sum(aer_lw_abs(:ncol, :, :), dim=3) + aodabslw(:ncol, :) = aer_lw_abs(:ncol, :, idx_lw_diag) + + !----------------------------------------------------------------- + ! Copy per-species fields for fillvalue modification + !----------------------------------------------------------------- + dustaod_l(:ncol) = dustaod(:ncol) + sulfaod_l(:ncol) = sulfaod(:ncol) + bcaod_l(:ncol) = bcaod(:ncol) + pomaod_l(:ncol) = pomaod(:ncol) + soaaod_l(:ncol) = soaaod(:ncol) + ssltaod_l(:ncol) = ssltaod(:ncol) + aodabsbc_l(:ncol) = aodabsbc(:ncol) + burdendust_l(:ncol) = burdendust(:ncol) + burdenso4_l(:ncol) = burdenso4(:ncol) + burdenbc_l(:ncol) = burdenbc(:ncol) + burdenpom_l(:ncol) = burdenpom(:ncol) + burdensoa_l(:ncol) = burdensoa(:ncol) + burdenseasalt_l(:ncol) = burdenseasalt(:ncol) + + !----------------------------------------------------------------- + ! Output day+night (dn) variants first (all columns valid) + !----------------------------------------------------------------- + + ! SW column + call history_out_field('AODVISdn', aodvis_l) + call history_out_field('AODUVdn', aoduv) + call history_out_field('AODNIRdn', aodnir) + call history_out_field('AODABSdn', aodabs) + call history_out_field('AODTOTdn', aodtot) + call history_out_field('SSAVISdn', ssavis_l) + + ! Stratospheric + call history_out_field('AODVISstdn', aodvisst) + call history_out_field('AODUVstdn', aoduvst) + call history_out_field('AODNIRstdn', aodnirst) + + ! SW per-level + call history_out_field('EXTINCTdn', extinct) + call history_out_field('EXTINCTUVdn', extinctuv) + call history_out_field('EXTINCTNIRdn', extinctnir) + call history_out_field('ABSORBdn', absorb) + + ! Per-species column (dn) + call history_out_field('AODDUSTdn', dustaod_l) + call history_out_field('AODSO4dn', sulfaod_l) + call history_out_field('AODBCdn', bcaod_l) + call history_out_field('AODPOMdn', pomaod_l) + call history_out_field('AODSOAdn', soaaod_l) + call history_out_field('AODSSdn', ssltaod_l) + call history_out_field('AODABSBCdn', aodabsbc_l) + + call history_out_field('BURDENDUSTdn', burdendust_l) + call history_out_field('BURDENSO4dn', burdenso4_l) + call history_out_field('BURDENBCdn', burdenbc_l) + call history_out_field('BURDENPOMdn', burdenpom_l) + call history_out_field('BURDENSOAdn', burdensoa_l) + call history_out_field('BURDENSEASALTdn', burdenseasalt_l) + + ! LW per-level (no day/night distinction) + call history_out_field('TOTABSLW', totabslw) + call history_out_field('AODABSLW', aodabslw) + + !----------------------------------------------------------------- + ! Apply fillvalue to night columns for day-only variants + !----------------------------------------------------------------- + do i = 1, nnite + aodvis_l(idxnite(i)) = fillvalue + aoduv(idxnite(i)) = fillvalue + aodnir(idxnite(i)) = fillvalue + aodabs(idxnite(i)) = fillvalue + aodtot(idxnite(i)) = fillvalue + ssavis_l(idxnite(i)) = fillvalue + aodvisst(idxnite(i)) = fillvalue + aoduvst(idxnite(i)) = fillvalue + aodnirst(idxnite(i)) = fillvalue + + extinct(idxnite(i), :) = fillvalue + extinctuv(idxnite(i), :) = fillvalue + extinctnir(idxnite(i), :) = fillvalue + absorb(idxnite(i), :) = fillvalue + + dustaod_l(idxnite(i)) = fillvalue + sulfaod_l(idxnite(i)) = fillvalue + bcaod_l(idxnite(i)) = fillvalue + pomaod_l(idxnite(i)) = fillvalue + soaaod_l(idxnite(i)) = fillvalue + ssltaod_l(idxnite(i)) = fillvalue + aodabsbc_l(idxnite(i)) = fillvalue + + burdendust_l(idxnite(i)) = fillvalue + burdenso4_l(idxnite(i)) = fillvalue + burdenbc_l(idxnite(i)) = fillvalue + burdenpom_l(idxnite(i)) = fillvalue + burdensoa_l(idxnite(i)) = fillvalue + burdenseasalt_l(idxnite(i)) = fillvalue + end do + + !----------------------------------------------------------------- + ! Output day-only variants (night columns have fillvalue) + !----------------------------------------------------------------- + + ! SW column + call history_out_field('AODVIS', aodvis_l) + call history_out_field('AODUV', aoduv) + call history_out_field('AODNIR', aodnir) + call history_out_field('AODABS', aodabs) + call history_out_field('AODTOT', aodtot) + call history_out_field('SSAVIS', ssavis_l) + + ! Stratospheric + call history_out_field('AODVISst', aodvisst) + call history_out_field('AODUVst', aoduvst) + call history_out_field('AODNIRst', aodnirst) + + ! SW per-level + call history_out_field('EXTINCT', extinct) + call history_out_field('EXTINCTUV', extinctuv) + call history_out_field('EXTINCTNIR', extinctnir) + call history_out_field('ABSORB', absorb) + + ! Per-species column + call history_out_field('AODDUST', dustaod_l) + call history_out_field('AODSO4', sulfaod_l) + call history_out_field('AODBC', bcaod_l) + call history_out_field('AODPOM', pomaod_l) + call history_out_field('AODSOA', soaaod_l) + call history_out_field('AODSS', ssltaod_l) + call history_out_field('AODABSBC', aodabsbc_l) + + call history_out_field('BURDENDUST', burdendust_l) + call history_out_field('BURDENSO4', burdenso4_l) + call history_out_field('BURDENBC', burdenbc_l) + call history_out_field('BURDENPOM', burdenpom_l) + call history_out_field('BURDENSOA', burdensoa_l) + call history_out_field('BURDENSEASALT', burdenseasalt_l) + + ! AEROD_v and AODvstrt (aliases for AODVIS/AODVISst, day-only) + call history_out_field('AEROD_v', aodvis_l) + call history_out_field('AODvstrt', aodvisst) + + ! Per-aerosol visible OD (ODV_), day-only with fillvalue + do iaer = 1, min(num_bulk_aer, num_odv) + odv_tmp(:ncol) = odv_col_aod(:ncol, iaer) + do i = 1, nnite + odv_tmp(idxnite(i)) = fillvalue + end do + call history_out_field(trim(odv_names(iaer)), odv_tmp) + end do + + ! Output diagnostics for active climate and diagnostic lists + do ilist = 0, N_DIAG + if (active_calls(ilist)) then + call rad_aer_diag_out(ilist, constituents, pdeldry, ncol) + end if + end do + + end subroutine aerosol_optics_diagnostics_run + +end module aerosol_optics_diagnostics diff --git a/schemes/sima_diagnostics/aerosol_optics_diagnostics.meta b/schemes/sima_diagnostics/aerosol_optics_diagnostics.meta new file mode 100644 index 00000000..abc7622c --- /dev/null +++ b/schemes/sima_diagnostics/aerosol_optics_diagnostics.meta @@ -0,0 +1,269 @@ +[ccpp-table-properties] + name = aerosol_optics_diagnostics + type = scheme + +[ccpp-arg-table] + name = aerosol_optics_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = aerosol_optics_diagnostics_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ nswbands ] + standard_name = number_of_bands_for_shortwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ nlwbands ] + standard_name = number_of_bands_for_longwave_radiation + units = count + type = integer + dimensions = () + intent = in +[ N_DIAG ] + standard_name = number_of_radiative_aerosol_diagnostic_lists + units = count + type = integer + dimensions = () + intent = in +[ active_calls ] + standard_name = flag_for_active_radiative_aerosol_diagnostic_list + units = flag + type = logical + dimensions = (index_of_climate_radiative_aerosol_list:number_of_radiative_aerosol_diagnostic_lists) + intent = in +[ constituents ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = in +[ aer_tau ] + standard_name = shortwave_aerosol_extinction_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = in +[ aer_tau_w ] + standard_name = shortwave_aerosol_single_scattering_albedo_times_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_shortwave_radiation) + intent = in +[ aer_lw_abs ] + standard_name = longwave_aerosol_absorption_optical_depth + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bands_for_longwave_radiation) + intent = in +[ idx_sw_diag ] + standard_name = index_of_shortwave_band + units = index + type = integer + dimensions = () + intent = in +[ idx_uv_diag ] + standard_name = index_of_UV_band + units = index + type = integer + dimensions = () + intent = in +[ idx_nir_diag ] + standard_name = index_of_near_IR_band + units = index + type = integer + dimensions = () + intent = in +[ idx_lw_diag ] + standard_name = index_of_longwave_band + units = index + type = integer + dimensions = () + intent = in +[ pdeldry ] + standard_name = air_pressure_thickness_of_dry_air + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rga ] + standard_name = reciprocal_of_gravitational_acceleration + units = s2 m-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ nnite ] + standard_name = nighttime_columns_dimension + units = count + type = integer + dimensions = () + intent = in +[ idxnite ] + standard_name = nighttime_columns + units = index + type = integer + dimensions = (horizontal_loop_extent) + intent = in +[ troplev ] + standard_name = tropopause_vertical_layer_index_from_chemical_method + units = index + type = integer + dimensions = (horizontal_loop_extent) + intent = in +[ dustaod ] + standard_name = aerosol_optical_depth_from_dust + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ sulfaod ] + standard_name = aerosol_optical_depth_from_sulfate + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ bcaod ] + standard_name = aerosol_optical_depth_from_black_carbon + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ pomaod ] + standard_name = aerosol_optical_depth_from_primary_organic_matter + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ soaaod ] + standard_name = aerosol_optical_depth_from_secondary_organic_aerosol + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ssltaod ] + standard_name = aerosol_optical_depth_from_seasalt + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ aodabsbc ] + standard_name = absorption_aerosol_optical_depth_from_black_carbon + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdendust ] + standard_name = column_burden_from_dust + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdenso4 ] + standard_name = column_burden_from_sulfate + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdenbc ] + standard_name = column_burden_from_black_carbon + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdenpom ] + standard_name = column_burden_from_primary_organic_matter + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdensoa ] + standard_name = column_burden_from_secondary_organic_aerosol + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ burdenseasalt ] + standard_name = column_burden_from_seasalt + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ssavis ] + standard_name = aerosol_single_scatter_albedo_at_vis + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ aodvis ] + standard_name = aerosol_optical_depth_at_vis + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ num_bulk_aer ] + standard_name = number_of_bulk_aerosol_constituents_in_climate_radiative_aerosol_list + units = count + type = integer + dimensions = () + intent = in +[ odv_col_aod ] + standard_name = column_visible_optical_depth_per_bulk_aerosol + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, number_of_bulk_aerosol_constituents_in_climate_radiative_aerosol_list) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/ndrop_bam_diagnostics.F90 b/schemes/sima_diagnostics/ndrop_bam_diagnostics.F90 new file mode 100644 index 00000000..669348b1 --- /dev/null +++ b/schemes/sima_diagnostics/ndrop_bam_diagnostics.F90 @@ -0,0 +1,87 @@ +! Diagnostics for ndrop_bam_ccpp scheme (BAM droplet activation) +! +! Requires ndrop_bam_ccpp_init to have run first (populates aername). +! This ordering is guaranteed by SDF placement: diagnostic schemes +! follow the physics schemes they diagnose. +module ndrop_bam_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: ndrop_bam_diagnostics_init + public :: ndrop_bam_diagnostics_run + + integer :: num_aer = 0 + +contains + + !> \section arg_table_ndrop_bam_diagnostics_init Argument Table + !! \htmlinclude ndrop_bam_diagnostics_init.html + subroutine ndrop_bam_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + use ndrop_bam, only: aername + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: l + + errmsg = '' + errflg = 0 + + ! CCN concentrations at fixed supersaturation levels + call history_add_field('CCN1', 'cloud_condensation_nuclei_number_concentration_at_S_0.02pct', 'lev', 'avg', 'cm-3') + call history_add_field('CCN2', 'cloud_condensation_nuclei_number_concentration_at_S_0.05pct', 'lev', 'avg', 'cm-3') + call history_add_field('CCN3', 'cloud_condensation_nuclei_number_concentration_at_S_0.1pct', 'lev', 'avg', 'cm-3') + call history_add_field('CCN4', 'cloud_condensation_nuclei_number_concentration_at_S_0.2pct', 'lev', 'avg', 'cm-3') + call history_add_field('CCN5', 'cloud_condensation_nuclei_number_concentration_at_S_0.5pct', 'lev', 'avg', 'cm-3') + call history_add_field('CCN6', 'cloud_condensation_nuclei_number_concentration_at_S_1.0pct', 'lev', 'avg', 'cm-3') + + ! Per-aerosol number concentration diagnostics (dynamic names from ndrop_bam) + if (allocated(aername)) then + num_aer = size(aername) + do l = 1, num_aer + call history_add_field(trim(aername(l))//'_m3', & + 'aerosol_number_concentration_of_'//trim(aername(l)), 'lev', 'avg', 'm-3') + end do + end if + + end subroutine ndrop_bam_diagnostics_init + + !> \section arg_table_ndrop_bam_diagnostics_run Argument Table + !! \htmlinclude ndrop_bam_diagnostics_run.html + subroutine ndrop_bam_diagnostics_run( & + psat, naer_all, & + ccn, naer2_diag, & + errmsg, errflg) + + use cam_history, only: history_out_field + use ndrop_bam, only: aername, ccn_name + + integer, intent(in) :: psat + integer, intent(in) :: naer_all + real(kind_phys), intent(in) :: ccn(:,:,:) + real(kind_phys), intent(in) :: naer2_diag(:,:,:) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: l + + errmsg = '' + errflg = 0 + + ! CCN concentrations (sliced from 3D array) + do l = 1, psat + call history_out_field(ccn_name(l), ccn(:,:,l)) + end do + + ! Per-aerosol number concentration diagnostics + do l = 1, min(naer_all, num_aer) + call history_out_field(trim(aername(l))//'_m3', naer2_diag(:,:,l)) + end do + + end subroutine ndrop_bam_diagnostics_run + +end module ndrop_bam_diagnostics diff --git a/schemes/sima_diagnostics/ndrop_bam_diagnostics.meta b/schemes/sima_diagnostics/ndrop_bam_diagnostics.meta new file mode 100644 index 00000000..dc69d822 --- /dev/null +++ b/schemes/sima_diagnostics/ndrop_bam_diagnostics.meta @@ -0,0 +1,60 @@ +[ccpp-table-properties] + name = ndrop_bam_diagnostics + type = scheme + dependencies = ../microp_aero/ndrop_bam.F90 + +[ccpp-arg-table] + name = ndrop_bam_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = ndrop_bam_diagnostics_run + type = scheme +[ psat ] + standard_name = number_of_ccn_supersaturation_levels + units = count + type = integer + dimensions = () + intent = in +[ naer_all ] + standard_name = number_of_bulk_aerosol_types + units = count + type = integer + dimensions = () + intent = in +[ ccn ] + standard_name = cloud_condensation_nuclei_number_concentration_at_supersaturation + units = cm-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccn_supersaturation_levels) + intent = in +[ naer2_diag ] + standard_name = bulk_aerosol_number_concentration_diagnostic + units = m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_bulk_aerosol_types) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/nucleate_ice_diagnostics.F90 b/schemes/sima_diagnostics/nucleate_ice_diagnostics.F90 new file mode 100644 index 00000000..115c4f61 --- /dev/null +++ b/schemes/sima_diagnostics/nucleate_ice_diagnostics.F90 @@ -0,0 +1,119 @@ +! Diagnostics for nucleate_ice_ccpp scheme +module nucleate_ice_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: nucleate_ice_diagnostics_init + public :: nucleate_ice_diagnostics_run + +contains + + !> \section arg_table_nucleate_ice_diagnostics_init Argument Table + !! \htmlinclude nucleate_ice_diagnostics_init.html + subroutine nucleate_ice_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + ! Ice nucleation number tendencies by mechanism + call history_add_field('NIHFTEN', 'tendency_of_ice_nucleation_number_concentration_from_homogeneous_freezing', & + 'lev', 'avg', 'm-3 s-1') + call history_add_field('NIIMMTEN', 'tendency_of_ice_nucleation_number_concentration_from_immersion_freezing', & + 'lev', 'avg', 'm-3 s-1') + call history_add_field('NIDEPTEN', 'tendency_of_ice_nucleation_number_concentration_from_deposition_nucleation', & + 'lev', 'avg', 'm-3 s-1') + call history_add_field('NIMEYTEN', 'tendency_of_ice_nucleation_number_concentration_from_meyers_nucleation', & + 'lev', 'avg', 'm-3 s-1') + + ! Regime, subgrid, and tropopause diagnostics + call history_add_field('NIREGM', 'ice_nucleation_regime_temperature_threshold', 'lev', 'avg', 'degC') + call history_add_field('NISUBGRID', 'subgrid_saturation_scaling_factor_diagnostic_for_ice_nucleation', 'lev', 'avg', '1') + call history_add_field('NITROP_PD', 'vertical_layer_indicator_of_chemical_tropopause', 'lev', 'avg', '1') + + ! Pre-existing ice diagnostics + call history_add_field('FHOM', 'fraction_of_cirrus_with_homogeneous_freezing', 'lev', 'avg', 'fraction') + call history_add_field('WICE', 'vertical_velocity_reduction_from_preexisting_ice_for_ice_nucleation', 'lev', 'avg', 'm s-1') + call history_add_field('WEFF', 'effective_vertical_velocity_for_ice_nucleation', 'lev', 'avg', 'm s-1') + + ! Aerosol input number concentration tendencies to ice nucleation + call history_add_field('INnso4TEN', 'tendency_of_sulfate_number_concentration_input_to_ice_nucleation', 'lev', 'avg', 'm-3 s-1') + call history_add_field('INnbcTEN', 'tendency_of_black_carbon_number_concentration_input_to_ice_nucleation', 'lev', 'avg', 'm-3 s-1') + call history_add_field('INndustTEN', 'tendency_of_dust_number_concentration_input_to_ice_nucleation', 'lev', 'avg', 'm-3 s-1') + call history_add_field('INondustTEN', 'tendency_of_dust_number_concentration_output_from_ice_nucleation', 'lev', 'avg', 'm-3 s-1') + call history_add_field('INhetTEN', 'tendency_of_heterogeneous_ice_nuclei_number_concentration', 'lev', 'avg', 'm-3 s-1') + call history_add_field('INhomTEN', 'tendency_of_homogeneous_ice_nuclei_number_concentration', 'lev', 'avg', 'm-3 s-1') + + ! Ice nucleation frequency diagnostics + call history_add_field('INFrehom', 'frequency_of_homogeneous_ice_nucleation', 'lev', 'avg', '1') + call history_add_field('INFreIN', 'frequency_of_ice_nucleation', 'lev', 'avg', '1') + + end subroutine nucleate_ice_diagnostics_init + + !> \section arg_table_nucleate_ice_diagnostics_run Argument Table + !! \htmlinclude nucleate_ice_diagnostics_run.html + subroutine nucleate_ice_diagnostics_run( & + nihf, niimm, nidep, nimey, & + regm, subgrid_diag, trop_pd, & + fhom, wice, weff, & + INnso4, INnbc, INndust, INondust, & + INhet, INhom, INFrehom, INFreIN, & + errmsg, errflg) + + use cam_history, only: history_out_field + + real(kind_phys), intent(in) :: nihf(:,:) + real(kind_phys), intent(in) :: niimm(:,:) + real(kind_phys), intent(in) :: nidep(:,:) + real(kind_phys), intent(in) :: nimey(:,:) + real(kind_phys), intent(in) :: regm(:,:) + real(kind_phys), intent(in) :: subgrid_diag(:,:) + real(kind_phys), intent(in) :: trop_pd(:,:) + real(kind_phys), intent(in) :: fhom(:,:) + real(kind_phys), intent(in) :: wice(:,:) + real(kind_phys), intent(in) :: weff(:,:) + real(kind_phys), intent(in) :: INnso4(:,:) + real(kind_phys), intent(in) :: INnbc(:,:) + real(kind_phys), intent(in) :: INndust(:,:) + real(kind_phys), intent(in) :: INondust(:,:) + real(kind_phys), intent(in) :: INhet(:,:) + real(kind_phys), intent(in) :: INhom(:,:) + real(kind_phys), intent(in) :: INFrehom(:,:) + real(kind_phys), intent(in) :: INFreIN(:,:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + call history_out_field('NIHFTEN', nihf) + call history_out_field('NIIMMTEN', niimm) + call history_out_field('NIDEPTEN', nidep) + call history_out_field('NIMEYTEN', nimey) + + call history_out_field('NIREGM', regm) + call history_out_field('NISUBGRID', subgrid_diag) + call history_out_field('NITROP_PD', trop_pd) + + call history_out_field('FHOM', fhom) + call history_out_field('WICE', wice) + call history_out_field('WEFF', weff) + + call history_out_field('INnso4TEN', INnso4) + call history_out_field('INnbcTEN', INnbc) + call history_out_field('INndustTEN', INndust) + call history_out_field('INondustTEN', INondust) + call history_out_field('INhetTEN', INhet) + call history_out_field('INhomTEN', INhom) + + call history_out_field('INFrehom', INFrehom) + call history_out_field('INFreIN', INFreIN) + + end subroutine nucleate_ice_diagnostics_run + +end module nucleate_ice_diagnostics diff --git a/schemes/sima_diagnostics/nucleate_ice_diagnostics.meta b/schemes/sima_diagnostics/nucleate_ice_diagnostics.meta new file mode 100644 index 00000000..a99ba135 --- /dev/null +++ b/schemes/sima_diagnostics/nucleate_ice_diagnostics.meta @@ -0,0 +1,143 @@ +[ccpp-table-properties] + name = nucleate_ice_diagnostics + type = scheme + +[ccpp-arg-table] + name = nucleate_ice_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = nucleate_ice_diagnostics_run + type = scheme +[ nihf ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_homogeneous_freezing + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ niimm ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_immersion_freezing + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ nidep ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_deposition_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ nimey ] + standard_name = tendency_of_ice_nucleation_number_concentration_from_meyers_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ regm ] + standard_name = ice_nucleation_regime_temperature_threshold + units = degC + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ subgrid_diag ] + standard_name = subgrid_saturation_scaling_factor_diagnostic_for_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ trop_pd ] + standard_name = vertical_layer_indicator_of_chemical_tropopause + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ fhom ] + standard_name = fraction_of_cirrus_with_homogeneous_freezing + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ wice ] + standard_name = vertical_velocity_reduction_from_preexisting_ice_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ weff ] + standard_name = effective_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INnso4 ] + standard_name = tendency_of_sulfate_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INnbc ] + standard_name = tendency_of_black_carbon_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INndust ] + standard_name = tendency_of_dust_number_concentration_input_to_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INondust ] + standard_name = tendency_of_dust_number_concentration_output_from_ice_nucleation + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INhet ] + standard_name = tendency_of_heterogeneous_ice_nuclei_number_concentration + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INhom ] + standard_name = tendency_of_homogeneous_ice_nuclei_number_concentration + units = m-3 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INFrehom ] + standard_name = frequency_of_homogeneous_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ INFreIN ] + standard_name = frequency_of_ice_nucleation + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.F90 b/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.F90 new file mode 100644 index 00000000..8432e1fe --- /dev/null +++ b/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.F90 @@ -0,0 +1,50 @@ +! Diagnostics for subgrid vertical velocity (WSUB, WSUBI) +module scale_subgrid_vertical_velocity_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: scale_subgrid_vertical_velocity_diagnostics_init + public :: scale_subgrid_vertical_velocity_diagnostics_run + +contains + + !> \section arg_table_scale_subgrid_vertical_velocity_diagnostics_init Argument Table + !! \htmlinclude scale_subgrid_vertical_velocity_diagnostics_init.html + subroutine scale_subgrid_vertical_velocity_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + call history_add_field('WSUB', 'subgrid_scale_vertical_velocity_for_droplet_nucleation', 'lev', 'avg', 'm s-1') + call history_add_field('WSUBI', 'subgrid_scale_vertical_velocity_for_ice_nucleation', 'lev', 'avg', 'm s-1') + + end subroutine scale_subgrid_vertical_velocity_diagnostics_init + + !> \section arg_table_scale_subgrid_vertical_velocity_diagnostics_run Argument Table + !! \htmlinclude scale_subgrid_vertical_velocity_diagnostics_run.html + subroutine scale_subgrid_vertical_velocity_diagnostics_run( & + wsub, wsubi, & + errmsg, errflg) + + use cam_history, only: history_out_field + + real(kind_phys), intent(in) :: wsub(:,:) + real(kind_phys), intent(in) :: wsubi(:,:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + call history_out_field('WSUB', wsub) + call history_out_field('WSUBI', wsubi) + + end subroutine scale_subgrid_vertical_velocity_diagnostics_run + +end module scale_subgrid_vertical_velocity_diagnostics diff --git a/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.meta b/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.meta new file mode 100644 index 00000000..ce0eec88 --- /dev/null +++ b/schemes/sima_diagnostics/scale_subgrid_vertical_velocity_diagnostics.meta @@ -0,0 +1,47 @@ +[ccpp-table-properties] + name = scale_subgrid_vertical_velocity_diagnostics + type = scheme + +[ccpp-arg-table] + name = scale_subgrid_vertical_velocity_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = scale_subgrid_vertical_velocity_diagnostics_run + type = scheme +[ wsub ] + standard_name = subgrid_scale_vertical_velocity_for_droplet_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ wsubi ] + standard_name = subgrid_scale_vertical_velocity_for_ice_nucleation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/vertical_diffusion/diffusion_solver.F90 b/schemes/vertical_diffusion/diffusion_solver.F90 index af3266bd..d4c5d0b9 100644 --- a/schemes/vertical_diffusion/diffusion_solver.F90 +++ b/schemes/vertical_diffusion/diffusion_solver.F90 @@ -847,8 +847,11 @@ subroutine vertical_diffusion_diffuse_tracers_init( & index(cnst_name, '_a2') > 0 .or. & ! MAM modal species. index(cnst_name, '_a3') > 0 .or. & ! MAM modal species. index(cnst_name, '_a4') > 0 .or. & ! MAM modal species. - index(cnst_name, 'NUMLIQ') > 0) then ! micro_pumas_cam if modal/trop_strat carma - ! is a modal aerosol species undergoing ndrop activation mixing. + index(cnst_name, 'NUMLIQ') > 0 .or. & ! micro_pumas_cam if modal/trop_strat carma + index(cnst_name, 'mass_number_concentration_of_cloud_liquid_wrt_moist_air_and_condensed_water') > 0) then + ! micro_pumas_cam (including MG1) if modal/trop_strat carma + + ! is a modal aerosol species (or numliq) undergoing ndrop activation mixing. ! in this case, vertical diffusion does not diffuse it. do_diffusion_const(m) = .false. end if diff --git a/suites/suite_cam4.xml b/suites/suite_cam4.xml index 41ac16cd..d6b09782 100644 --- a/suites/suite_cam4.xml +++ b/suites/suite_cam4.xml @@ -9,6 +9,8 @@ Microphysics RK Radiation RRTMGP Chemistry None (not implemented) + Aerosol Model BAM (bulk) + Aerosol Source Prescribed Vertical Diffusion HB Gravity Wave Drag Orographic --> @@ -20,6 +22,16 @@ zm_conv_options + + prescribed_ozone + prescribed_aerosols + prescribed_aerosol_deposition_flux + tropopause_find + prescribed_volcanic_aerosol + + + rebin_seasalt + check_energy_gmean @@ -215,8 +227,17 @@ sima_state_diagnostics - + rrtmgp_pre + + + aerosol_optics + aerosol_optics_diagnostics + + rrtmgp_cloud_optics_setup tropopause_find rrtmgp_variables diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index 76fb3e13..e1462fba 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -1,15 +1,42 @@ + - + to_be_ccppized_temporary + prescribe_radiative_gas_concentrations set_cloud_fraction_top + + prescribed_ozone + prescribed_aerosols + prescribed_aerosol_deposition_flux + tropopause_find + prescribed_volcanic_aerosol + + + rebin_seasalt + check_energy_gmean @@ -88,10 +115,167 @@ + + + + + + + + + + + + + tropopause_find + + + compute_subgrid_vertical_velocity_clubb + scale_subgrid_vertical_velocity + scale_subgrid_vertical_velocity_diagnostics + + + dust_default_radii + + + nucleate_ice_ccpp + nucleate_ice_diagnostics + + + ndrop_bam_ccpp + ndrop_bam_diagnostics + + + hetfrz_classnuc_stub + + + + + + + + + sima_state_diagnostics - + + rrtmgp_pre + + + aerosol_optics + aerosol_optics_diagnostics + + + rrtmgp_cloud_optics_setup + tropopause_find + rrtmgp_variables + rrtmgp_inputs + rrtmgp_sw_cloud_optics + rrtmgp_sw_mcica_subcol_gen + rrtmgp_cloud_diagnostics + + + rrtmgp_constituents + rrtmgp_sw_gas_optics_pre + rrtmgp_sw_gas_optics + solar_irradiance_data + rrtmgp_sw_solar_var + rrtmgp_sw_aerosols + rrtmgp_sw_rte + rrtmgp_sw_calculate_fluxes + rrtmgp_sw_calculate_heating_rate + rrtmgp_sw_diagnostics + rrtmgp_subcycle + + rrtmgp_lw_cloud_optics + rrtmgp_lw_mcica_subcol_gen + + + rrtmgp_constituents + rrtmgp_lw_gas_optics_pre + rrtmgp_lw_gas_optics + rrtmgp_lw_aerosols + rrtmgp_lw_rte + rrtmgp_lw_calculate_fluxes + rrtmgp_lw_calculate_heating_rate + rrtmgp_lw_diagnostics + rrtmgp_subcycle + + + rrtmgp_inputs_setup + rrtmgp_sw_solar_var_setup + rrtmgp_dry_static_energy_tendency + calculate_net_heating + rrtmgp_post + rrtmgp_diagnostics + + apply_heating_rate + geopotential_temp + + + + + @@ -99,9 +283,6 @@ tropopause_find tropopause_diagnostics - - - @@ -135,6 +316,9 @@ convert_dry_constituent_tendencies_to_dry_air_basis + + gravity_wave_drag_common_diagnostics + apply_tendency_of_eastward_wind apply_tendency_of_northward_wind @@ -144,6 +328,14 @@ geopotential_temp update_dry_static_energy + + check_energy_scaling + check_energy_chng + + + check_energy_save_teout diff --git a/test/test_schemes/initialize_pbuf_aerosols.F90 b/test/test_schemes/initialize_pbuf_aerosols.F90 new file mode 100644 index 00000000..cc073a6a --- /dev/null +++ b/test/test_schemes/initialize_pbuf_aerosols.F90 @@ -0,0 +1,136 @@ +! Initialize prescribed aerosols from file. +! This is the "pbuf" equivalent of initialize_constituents: +! when aerosols are prescribed (in a "-chem none" configuration in CAM) +! they are not constituents, but pbuf_ fields instead (rad_climate N: specifier) +! +! In CAM-SIMA, the "provider" of chemistry/aerosol needs to register these species, +! thus a scheme needs to register these constituents based on pbuf_ +! entries in the ncdata input file. +! +! This scheme should ONLY be used for snapshot runs. +! For non-snapshot runs, use the prescribed_aerosols and associated schemes. +module initialize_pbuf_aerosols + + implicit none + private + + public :: initialize_pbuf_aerosols_register + + integer, parameter :: num_bam_species = 13 + character(len=5), parameter :: bam_species(num_bam_species) = & + (/'sulf ', 'bcar1', 'bcar2', 'ocar1', 'ocar2', & + 'sslt1', 'sslt2', 'sslt3', 'sslt4', & + 'dust1', 'dust2', 'dust3', 'dust4'/) + + integer, parameter :: num_mam_species = 30 + character(len=6), parameter :: mam_species(num_mam_species) = & + (/'num_a1', 'bc_a1 ', 'dst_a1', 'ncl_a1', 'pom_a1', 'so4_a1', 'soa_a1', & + 'num_a2', 'ncl_a2', 'so4_a2', 'soa_a2', & + 'num_a3', 'dst_a3', 'ncl_a3', 'so4_a3', & + 'num_c1', 'bc_c1 ', 'dst_c1', 'ncl_c1', 'pom_c1', 'so4_c1', 'soa_c1', & + 'num_c2', 'ncl_c2', 'so4_c2', 'soa_c2', & + 'num_c3', 'dst_c3', 'ncl_c3', 'so4_c3'/) + +contains + +!> \section arg_table_initialize_pbuf_aerosols_register Argument Table +!! \htmlinclude initialize_pbuf_aerosols_register.html + subroutine initialize_pbuf_aerosols_register(constituents, errmsg, errflg) + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use cam_initfiles, only: initial_file_get_id + use pio, only: file_desc_t, var_desc_t, & + pio_inq_varid, pio_seterrorhandling, & + PIO_BCAST_ERROR, PIO_NOERR + use ccpp_kinds, only: kind_phys + use ccpp_chem_utils, only: chem_constituent_qmin + + ! Dummy variables + type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituents(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + type(file_desc_t), pointer :: ncdata + type(var_desc_t) :: varid + integer :: ierr + integer :: old_err_method + integer :: i, nspecies + logical :: is_bam, is_mam + character(len=6), allocatable :: species_list(:) + character(len=7) :: unit_name + + character(len=*), parameter :: subname = 'initialize_pbuf_aerosols_register' + + errflg = 0 + errmsg = '' + + ncdata => initial_file_get_id() + + ! Suppress PIO fatal errors for variable probing + call pio_seterrorhandling(ncdata, PIO_BCAST_ERROR, old_err_method) + + ! Detect aerosol mode by probing for indicator variables + is_mam = .false. + is_bam = .false. + + ierr = pio_inq_varid(ncdata, 'pbuf_num_a1', varid) + if (ierr == PIO_NOERR) is_mam = .true. + + ierr = pio_inq_varid(ncdata, 'pbuf_sulf', varid) + if (ierr == PIO_NOERR) is_bam = .true. + + ! Restore PIO error handling + call pio_seterrorhandling(ncdata, old_err_method) + + if (.not. is_bam .and. .not. is_mam) then + errflg = 1 + errmsg = subname//': No BAM (pbuf_sulf) or MAM (pbuf_num_a1) aerosol fields found on initial file. This scheme should NOT be used for non-snapshot runs or physics tests that do not use aerosol' + return + end if + + ! Build species list based on detected mode (MAM takes priority) + if (is_mam) then + nspecies = num_mam_species + allocate(species_list(nspecies)) + species_list(:) = mam_species(:) + else + nspecies = num_bam_species + allocate(species_list(nspecies)) + species_list(:) = bam_species(:) + end if + + ! Register constituents + allocate(constituents(nspecies), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + errflg = 1 + errmsg = subname // ': Failed to allocate constituents: ' // trim(errmsg) + return + end if + + do i = 1, nspecies + if (index(species_list(i), 'num_') == 1) then + unit_name = 'kg-1' + else + unit_name = 'kg kg-1' + end if + + call constituents(i)%instantiate( & + std_name = trim(species_list(i)), & + diag_name = trim(species_list(i)), & + long_name = 'prescribed aerosol ' // trim(species_list(i)), & + units = trim(unit_name), & + vertical_dim = 'vertical_layer_dimension', & + min_value = chem_constituent_qmin(trim(species_list(i))), & + advected = .false., & + water_species = .false., & + mixing_ratio_type = 'dry', & + errcode = errflg, & + errmsg = errmsg) + if (errflg /= 0) return + end do + + deallocate(species_list) + + end subroutine initialize_pbuf_aerosols_register + +end module initialize_pbuf_aerosols diff --git a/test/test_schemes/initialize_pbuf_aerosols.meta b/test/test_schemes/initialize_pbuf_aerosols.meta new file mode 100644 index 00000000..d0ba2fa5 --- /dev/null +++ b/test/test_schemes/initialize_pbuf_aerosols.meta @@ -0,0 +1,27 @@ +[ccpp-table-properties] + name = initialize_pbuf_aerosols + type = scheme +[ccpp-arg-table] + name = initialize_pbuf_aerosols_register + type = scheme +[ constituents ] + standard_name = dynamic_constituents_for_initialize_pbuf_aerosols + units = none + dimensions = (:) + allocatable = True + type = ccpp_constituent_properties_t + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out diff --git a/test/test_suites/suite_microp_aero_bam.xml b/test/test_suites/suite_microp_aero_bam.xml new file mode 100644 index 00000000..8b8096ff --- /dev/null +++ b/test/test_suites/suite_microp_aero_bam.xml @@ -0,0 +1,74 @@ + + + + + + + + initialize_constituents + + + initialize_pbuf_aerosols + + + to_be_ccppized_temporary + + + rebin_seasalt + + + tropopause_find + + + compute_subgrid_vertical_velocity_clubb + scale_subgrid_vertical_velocity + scale_subgrid_vertical_velocity_diagnostics + + + dust_default_radii + + + nucleate_ice_ccpp + nucleate_ice_diagnostics + + + ndrop_bam_ccpp + ndrop_bam_diagnostics + + + hetfrz_classnuc_stub + + diff --git a/test/test_suites/suite_rrtmgp.xml b/test/test_suites/suite_rrtmgp.xml index edc64652..00d2a7cd 100644 --- a/test/test_suites/suite_rrtmgp.xml +++ b/test/test_suites/suite_rrtmgp.xml @@ -11,6 +11,7 @@ rrtmgp_sw_cloud_optics rrtmgp_sw_mcica_subcol_gen rrtmgp_cloud_diagnostics + aerosol_optics_zero_stub rrtmgp_constituents