diff --git a/Externals_CAM.cfg b/Externals_CAM.cfg index ce4ed313ce..73422bf06f 100644 --- a/Externals_CAM.cfg +++ b/Externals_CAM.cfg @@ -8,8 +8,8 @@ required = True [carma] local_path = src/physics/carma/base protocol = git -repo_url = https://github.com/ESCOMP/CARMA_base.git -tag = carma4_05 +repo_url = https://github.com/fvitt/CARMA_base.git +branch = carma4_05_concelem required = True [cosp2] diff --git a/src/physics/carma/cam/carma_intr.F90 b/src/physics/carma/cam/carma_intr.F90 index 1ef7c30543..9461c57322 100644 --- a/src/physics/carma/cam/carma_intr.F90 +++ b/src/physics/carma/cam/carma_intr.F90 @@ -70,6 +70,27 @@ module carma_intr public :: carma_restart_write public :: carma_restart_read + + ! Microphysics info from CAM state + ! + ! NOTE: These calls can be used in CAM when the CAM state is available, but the CARMASTATE + ! is not available. These will return the instantaneous values instead of relying on + ! pbuf fields that might be from the previous timestep. + public carma_get_bin + public carma_get_bin_cld + public carma_get_dry_radius + public carma_get_elem_for_group + public carma_get_group_by_name + public carma_get_kappa + public carma_get_number + public carma_get_number_cld + public carma_get_total_mmr + public carma_get_total_mmr_cld + public carma_get_wet_radius + public carma_get_bin_rmass + public carma_set_bin + + ! NOTE: This is required by physpkg.F90, since the carma_intr.F90 stub in physics/cam ! does not have access to carma_constant.F90, but needs to also provide a defintion ! for MAXCLDAERDIAG. Thus the definition of this variable needs to come from @@ -944,6 +965,7 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli real(r8) :: ar(pver) ! Area Ratio real(r8) :: vm(pver) ! Massweighted fall velocity (cm2) real(r8) :: jn(pver) ! nucleation (cm-3) + real(r8) :: totalmmr(pver) ! total particle mmr (kg/kg) real(r8) :: numberDensity(pver) ! number density (cm-3) real(r8) :: nucleationRate(pver) ! nucleation rate (cm-3 s-1) real(r8) :: extinctionCoefficient(pver) ! extinction coefficient (cm2) @@ -1324,7 +1346,8 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli do ibin = 1, NBIN call CARMASTATE_GetBin(cstate, ielem, ibin, newstate(:), rc, & - numberDensity=numberDensity, nucleationRate=nucleationRate, r_wet=r_wet, rhop_wet=rhop_wet, sedimentationflux=dd, vd=vd, vf=vf, dtpart=dtpart) + numberDensity=numberDensity, nucleationRate=nucleationRate, r_wet=r_wet, & + rhop_wet=rhop_wet, sedimentationflux=dd, vd=vd, vf=vf, dtpart=dtpart, totalmmr=totalmmr) if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_GetBin failed.') ! For prognostic groups, set the tendency from the corresponding constituents. @@ -1355,7 +1378,7 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli re3(:) = re3(:) + numberDensity(:) * ((r(ibin)*rrat(ibin))**3) ad(:) = ad(:) + numberDensity(:) * 4.0_r8 * PI * (r(ibin)**2) * 1.0e8_r8 md(:) = md(:) + numberDensity(:) * rmass(ibin) - mr(:) = mr(:) + newstate(:) + mr(:) = mr(:) + totalmmr(:) pa(:) = pa(:) + numberDensity(:) * PI * ((r(ibin) * rrat(ibin))**2) * arat(ibin) vm(:) = vm(:) + numberDensity(:) * rmass(ibin) * vf(2:) / 100._f @@ -1370,8 +1393,8 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli od(:) = od(:) + numberDensity(:) * extinctionCoefficient(:) * dz(:) * 100._r8 end if - bndiags(icol,:,ibin,ielem,BNDIAGS_VR) = bndiags(icol,:,ibin,ielem,BNDIAGS_VR) + newstate(:) - gpdiags(icol, :, igroup, GPDIAGS_VR) = gpdiags(icol, :, igroup, GPDIAGS_VR) + newstate(:) + bndiags(icol,:,ibin,ielem,BNDIAGS_VR) = bndiags(icol,:,ibin,ielem,BNDIAGS_VR) + totalmmr(:) + gpdiags(icol, :, igroup, GPDIAGS_VR) = gpdiags(icol, :, igroup, GPDIAGS_VR) + totalmmr(:) ! Particle temperatures from particle heating. if (carma_do_pheat) then @@ -3684,4 +3707,721 @@ subroutine CARMA_restart_read(File) end subroutine CARMA_restart_read + + + !! Get the mixing ratio for the specified element and bin. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_bin(state, ielem, ibin, mmr, rc) + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: ielem !! element index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: mmr(pcols,pver) !! mass mixing ratio (kg/kg) + integer, intent(out) :: rc !! return code + + integer :: ncol + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the group and bin ranges + if ((ielem < 1) .or. (ielem .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_bin:: ERROR - Invalid element id, ', ielem + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_bin:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Get the element from the physics state + mmr(:ncol, :) = state%q(:ncol, :, icnst4elem(ielem, ibin)) + + return + end subroutine + + !! Get the mixing ratio for the specified element and bin. + subroutine carma_get_bin_cld(pbuf, ielem, ibin, ncol, nlev, mmr, rc) + type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer + integer, intent(in) :: ielem !! element index + integer, intent(in) :: ibin !! bin index + integer, intent(in) :: ncol,nlev !! dimensions + real(r8), intent(out) :: mmr(:,:) !! mass mixing ratio (kg/kg) + integer, intent(out) :: rc !! return code + + real(r8), pointer :: mmr_ptr(:,:) + character(len=8) :: shortname ! short (CAM) name + character(len=16) :: c_name + integer :: idx + + ! default return code + rc = RC_OK + + call CARMAELEMENT_Get(carma, ielem, rc, shortname=shortname) + + write(c_name, '(A, I2.2)') trim(shortname), ibin + + idx = pbuf_get_index('CLD'//trim(c_name)) + call pbuf_get_field(pbuf, idx, mmr_ptr) + + mmr(:ncol,:nlev) = mmr_ptr(:ncol,:nlev) + + end subroutine carma_get_bin_cld + + !! Determine the dry radius and dry density for the particular bin. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc) + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: rdry(:,:) !! dry radius (m) + real(r8), intent(out) :: rhopdry(:,:) !! dry density (kg/m3) + integer, intent(out) :: rc !! return code + + real(r8) :: rhoelem(NBIN) ! element density (g/cm3) + real(r8) :: totvol(pcols,pver) ! total volume (m3/kg) + real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg) + real(r8) :: mmr(pcols, pver) ! mass mixing ratio (kg/kg) + real(r8) :: nmr(pcols, pver) ! number mixing ratio (#/kg) + integer :: nelems ! number of elements in group + integer :: ielems(NELEM) ! element indexes for group + integer :: ncol + integer :: i + integer :: ielem + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_dry_radius:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_dry_radius:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Iterate over all of the composition and determine the dry volume and dry radius. + call carma_get_elem_for_group(igroup, nelems, ielems, rc) + if (rc < 0) return + + totvol(:ncol, :) = 0._r8 + totmmr(:ncol, :) = 0._r8 + + do i = 1, nelems + ielem = ielems(i) + + call CARMAELEMENT_Get(carma, ielem, rc, rho=rhoelem) + if (rc < 0) return + + call carma_get_bin(state, ielem, ibin, mmr, rc) + if (rc < 0) return + + totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :) + totvol(:ncol, :) = totvol(:ncol, :) + mmr(:ncol, :) / (rhoelem(ibin) / 1.e3_r8 * 1.e6_r8) + end do + + ! Add checks for totvol = 0 and nmr = 0 + rhopdry(:ncol, :) = totmmr(:ncol, :) / totvol(:ncol, :) + + call carma_get_number(state, igroup, ibin, nmr, rc) + if (rc < 0) return + + rdry(:ncol, :) = ((3._r8 * totvol(:ncol, :) / nmr(:ncol, :)) / (4._r8 * PI)) ** (1._r8 / 3._r8) + + return + end subroutine carma_get_dry_radius + + + !! Get the number of elements and list of element ids for a group. This includes + !! the concentration elements and the core masses. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_elem_for_group(igroup, nelems, ielems, rc) + integer, intent(in) :: igroup !! group index + integer, intent(out) :: nelems !! number of elements in group + integer, intent(out) :: ielems(NELEM) !! indexes of elements in group + integer, intent(out) :: rc !! return code + + integer :: ienconc + integer :: ncore + integer :: icorelem(NELEM) + + ! default return code + rc = RC_OK + + ! Check the group range. + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_elem_for_group:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + call CARMAGROUP_Get(carma, igroup, rc, ienconc=ienconc, ncore=ncore, icorelem=icorelem) + + nelems = ncore + 1 + ielems(1) = ienconc + + if (ncore .gt. 0) then + ielems(2:ncore+1) = icorelem(1:ncore) + end if + + return + end subroutine + + + !! Get the CARMA group id a group name. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_group_by_name(shortname, igroup, rc) + character(len=*), intent(in) :: shortname !! the group short name + integer, intent(out) :: igroup !! group index + integer, intent(out) :: rc !! return code + + integer :: i + character(len=32) :: name + + ! default return code + rc = RC_OK + + igroup = -1 + + ! Check the short names of each group for one that matches + do i = 1, NGROUP + call CARMAGROUP_Get(carma, i, rc, shortname=name) + + if (trim(shortname) .eq. trim(name)) then + igroup = i + exit + end if + end do + + if (igroup .eq. -1) then + write(LUNOPRT, *) 'carma_get_group_by_name:: ERROR - group not found, ', shortname + rc = RC_ERROR + return + end if + + return + end subroutine + + + !! Get the CARMA group id and bin id from a compound name xxxxxxnn, where xxxxxx is the + !! name of the group and nn is the two digit bin number. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_group_and_bin_by_name(shortname, igroup, ibin, rc) + character(len=*), intent(out) :: shortname !! the group short name + integer, intent(out) :: igroup !! group index + integer, intent(out) :: ibin !! bin index + integer, intent(out) :: rc !! return code + + integer :: i + character(len=32) :: name + character(len=32) :: groupname + character(len=32) :: binname + + ! default return code + rc = RC_OK + + igroup = -1 + ibin = -1 + + if (len(shortname) <= 2) then + write(LUNOPRT, *) 'carma_get_group_and_bin_by_name:: ERROR - Illegal shortname, ' // shortname + rc = RC_ERROR + return + end if + + ! Check the short names of each group for one that matches + groupname = shortname(:len(shortname)-2) + binname = shortname(len(shortname)-2:) + + call carma_get_group_by_name(groupname, igroup, rc) + if (rc < 0) return + + read(binname, *) ibin + + return + end subroutine + + + !! Determine a mass weighted kappa for the entire particle. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_kappa(state, igroup, ibin, kappa, rc) + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: kappa(pcols,pver) !! kappa value for the entire particle + integer, intent(out) :: rc !! return code + + real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg) + real(r8) :: mmr(pcols,pver) ! element mmr (kg/kg) + real(r8) :: kappaelem ! element kappa + integer :: ncol + integer :: nelems + integer :: ielems(NELEM) + integer :: i + integer :: ielem + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_kappa:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (igroup .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_kappa:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Iterate over all of the composition and determine the total mass. + call carma_get_elem_for_group(igroup, nelems, ielems, rc) + if (rc < 0) return + + totmmr(:ncol, :) = 0._r8 + kappa(:ncol, :) = 0._r8 + + do i = 1, nelems + ielem = ielems(i) + + call carma_get_bin(state, ielem, ibin, mmr, rc) + if (rc < 0) return + + call CARMAELEMENT_Get(carma, ielem, rc, kappa=kappaelem) + + kappa(:ncol, :) = kappa(:ncol, :) + mmr(:ncol, :) * kappaelem + totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :) + end do + + ! Figure out the average kappa.q + where (totmmr .gt. 0._r8) + kappa = kappa / totmmr + end where + + return + end subroutine + + + !! Get the number mixing ratio for the group. This is the number of particles per + !! density of air. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_number(state, igroup, ibin, nmr, rc) + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: nmr(pcols,pver) !! number mixing ratio (#/kg) + integer, intent(out) :: rc !! return code + + real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g) + real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg) + integer :: ncol + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (igroup .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Get the mass in each bin + call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) + if (rc < 0) return + + ! Get the total mmr in the bin + call carma_get_total_mmr(state, igroup, ibin, totmmr, rc) + if (rc < 0) return + + ! Get the mmr is the total mass divided by rmass, but need to convert rmass + ! to kg. + nmr(:ncol, :) = totmmr(:ncol, :) / (rmass(ibin) / 1.e3_r8) + + return + end subroutine carma_get_number + + subroutine carma_get_number_cld(pbuf, igroup, ibin, ncol, nlev, nmr, rc) + + implicit none + + type(physics_buffer_desc),pointer :: pbuf(:) !! physics buffer + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + integer, intent(in) :: ncol,nlev !! dimensions + real(r8), intent(out) :: nmr(pcols,pver) !! number mixing ratio (#/kg) + integer, intent(out) :: rc !! return code + + real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g) + real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg) + + ! default return code + rc = RC_OK + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (igroup .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Get the mass in each bin + call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) + if (rc < 0) return + + ! Get the total mmr in the bin + call carma_get_total_mmr_cld(pbuf, igroup, ibin, ncol, nlev, totmmr, rc) + if (rc < 0) return + + ! Get the mmr is the total mass divided by rmass, but need to convert rmass + ! to kg. + nmr(:ncol, :) = totmmr(:ncol, :) / (rmass(ibin) / 1.e3_r8) + + return + end subroutine carma_get_number_cld + + + !! Get the mixing ratio for the group. This is the total of all the elements that + !! make up the group. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_total_mmr(state, igroup, ibin, totmmr, rc) + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: totmmr(pcols,pver) !! total mmr (kg/kg) + integer, intent(out) :: rc !! return code + + real(r8) :: mmr(pcols, pver) ! mmr (kg/kg) + integer :: i + integer :: nelems + integer :: ielems(NELEM) + integer :: ielem + integer :: ncol + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Iterate over all of the composition and determine the total mass. + call carma_get_elem_for_group(igroup, nelems, ielems, rc) + if (rc < 0) return + + totmmr(:ncol, :) = 0._r8 + + do i = 1, nelems + ielem = ielems(i) + + call carma_get_bin(state, ielem, ibin, mmr, rc) + if (rc < 0) return + + totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :) + end do + + return + end subroutine carma_get_total_mmr + + subroutine carma_get_total_mmr_cld(pbuf, igroup, ibin, ncol, nlev, totmmr, rc) + + type(physics_buffer_desc),pointer :: pbuf(:) !! physics buffer + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + integer, intent(in) :: ncol,nlev !! dimensions + real(r8), intent(out) :: totmmr(pcols,pver) !! total mmr (kg/kg) + integer, intent(out) :: rc !! return code + + real(r8) :: mmr(pcols, pver) ! mmr (kg/kg) + integer :: i + integer :: nelems + integer :: ielems(NELEM) + integer :: ielem + + ! default return code + rc = RC_OK + + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Iterate over all of the composition and determine the total mass. + call carma_get_elem_for_group(igroup, nelems, ielems, rc) + if (rc < 0) return + + totmmr(:ncol, :) = 0._r8 + + do i = 1, nelems + ielem = ielems(i) + + call carma_get_bin_cld(pbuf, ielem, ibin, ncol, nlev, mmr, rc) + if (rc < 0) return + + totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :) + end do + + end subroutine carma_get_total_mmr_cld + + !! Find the wet radius and wet density for the group and bin specified. + !! + !! NOTE: Groups can be configured with different methods to determine the wet + !! radius, so multiple methods need to be supported and code from rhopart and + !! wetr need to be included in this routine. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) + use wv_sat_methods, only : wv_sat_qsat_water + use wetr, only : getwetr + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8), intent(out) :: rwet(pcols,pver) !! wet radius (m) + real(r8), intent(out) :: rhopwet(pcols,pver) !! wet density (kg/m3) + integer, intent(inout) :: rc !! return code + + real(r8) :: rdry(pcols,pver) !! dry radius (m) + real(r8) :: rhopdry(pcols,pver) !! dry density (kg/m3) + real(r8) :: rhoa(pcols,pver) !! air density (kg/m3) + real(r8) :: kappa(pcols,pver) !! dry radius (m) + real(r8) :: es !! saturation vapor pressure + real(r8) :: qs !! saturation specific humidity + real(r8) :: relhum !!relative humidity + integer :: icol + integer :: iz + integer :: ncol + integer :: iq + integer :: irhswell + + ! default return code + rc = RC_OK + + ncol = state%ncol + + + ! Check the group and bin ranges + if ((igroup < 1) .or. (igroup .gt. NGROUP)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Get the constiuent index for water vapor (Q) + call cnst_get_ind("Q", iq) + + ! The wet radius can be configured differently for each group, so we + ! need to use getwetr to handle those differences. This requires repeating + ! some code that is in rhopart to use getwetr properly. There may be a + ! better way to do this, but for now we will duplicate the code. + call carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc) + if (rc < 0) return + + ! Calculate the dry air density at each level, using the ideal gas law. + rhoa(:ncol, :) = (state%pmid(:ncol, :) * 10._r8) / (R_AIR * state%t(:ncol, :)) / 1.e3_r8 * 1.e6_r8 + + call CARMAGROUP_Get(carma, igroup, rc, irhswell=irhswell) + if (rc < 0) return + + do icol = 1, ncol + do iz = 1, pver + + ! Get relative humidity and vapor pressure + call wv_sat_qsat_water(state%t(icol,iz), state%pmid(icol,iz), es, qs) + relhum = state%q(icol,iz,iq) / qs + + ! If humidity affects the particle, then determine the equilbirium + ! radius and density based upon the relative humidity. + ! + ! NOTE: getwetr is in cgs units, so some conversions are needed from the + ! mks values + if (irhswell == I_WTPCT_H2SO4) then + + ! Should r be rdry and is rhop aready known? + call getwetr(carma, igroup, relhum, rdry(icol, iz) * 1e4_r8, rwet(icol, iz), & + rhopdry(icol, iz) * 1.e3_r8 / 1.e6_r8, rhopwet(icol,iz), rc, & + h2o_mass=state%q(icol,iz,iq) * rhoa(icol, iz) * 1.e3_r8 / 1.e6_r8, & + h2o_vp=es * 10._r8 / 1.e4_r8, temp=state%t(icol,iz)) + if (rc < 0) return + + else if (irhswell == I_PETTERS) then + call carma_get_kappa(state, igroup, ibin, kappa, rc) + if (rc < 0) return + + call getwetr(carma, igroup, relhum, rdry(icol, iz) * 1e4_r8, rwet(icol, iz), & + rhopdry(icol, iz) * 1.e3_r8 / 1.e6_r8, rhopwet(icol,iz), rc, & + h2o_mass=state%q(icol,iz,iq) * rhoa(icol, iz) * 1.e3_r8 / 1.e6_r8, & + h2o_vp=es * 10._r8 / 1.e4_r8, temp=state%t(icol,iz), kappa=kappa(icol,iz)) + if (rc < 0) return + + else ! I_GERBER and I_FITZGERALD + call getwetr(carma, igroup, relhum, rdry(icol, iz) * 1e4_r8, rwet(icol, iz), & + rhopdry(icol, iz) * 1.e3_r8 / 1.e6_r8, rhopwet(icol,iz), rc) + if (rc < 0) return + end if + end do + end do + + ! Convert rwet and rhopwet to mks units + rwet(:ncol,:) = rwet(:ncol,:) / 1.e4_r8 + rhopwet(:ncol,:) = rhopwet(:ncol,:) / 1.e3_r8 * 1.e6_r8 + + return + end subroutine + + + !! Provides the tendency (in kg/kg/s) required to change the element and bin from + !! the current state to the desired mmr. + !! + !! NOTE: The caller needs to make sure that the lq flags are set in ptend for the + !! particular tracer. Perhaps we need a routine that will set lq to true for all + !! the fields that could be set by CARMA to be used by the caller of this routine. + !! + !! @author Chuck Bardeen + !! @version Aug 2023 + subroutine carma_set_bin(state, ielem, ibin, mmr, dt, ptend, rc) + + implicit none + + type(physics_state), intent(in) :: state !! physics state variables + integer, intent(in) :: ielem !! element index + integer, intent(in) :: ibin !! bin index + real(r8), intent(in) :: mmr(pcols,pver) !! mass mixing ratio (kg/kg) + integer :: dt !! timestep size (sec) + type(physics_ptend), intent(inout) :: ptend !! constituent tendencies + integer, intent(out) :: rc !! return code + + integer :: ncol + integer :: icnst + + ! default return code + rc = RC_OK + + ncol = state%ncol + + ! Check the element and bin ranges + if ((ielem < 1) .or. (ielem .gt. NELEM)) then + write(LUNOPRT, *) 'carma_set_bin:: ERROR - Invalid element id, ', ielem + rc = RC_ERROR + return + end if + + if ((ibin < 1) .or. (ibin .gt. NBIN)) then + write(LUNOPRT, *) 'carma_set_binr:: ERROR - Invalid bin id, ', ibin + rc = RC_ERROR + return + end if + + ! Determine the tendency needed to make state into mmr for this tracer. + icnst = icnst4elem(ielem, ibin) + ptend%q(:ncol, :, icnst) = (mmr(:ncol, :) - state%q(:ncol, :, icnst)) / dt + + return + end subroutine + + subroutine carma_get_bin_rmass(igroup, ibin, mass, rc) + + integer, intent(in) :: igroup !! group index + integer, intent(in) :: ibin !! bin index + real(r8),intent(out) :: mass + integer, intent(out) :: rc !! return code + + real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g) + + ! default return code + rc = RC_OK + + call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) + if (rc /= RC_OK) return + + mass = rmass(ibin) + + end subroutine carma_get_bin_rmass + end module carma_intr diff --git a/src/physics/carma/models/trop_strat_soa1/carma_model_mod.F90 b/src/physics/carma/models/trop_strat_soa1/carma_model_mod.F90 index e4d6e8910d..4d69d1279f 100644 --- a/src/physics/carma/models/trop_strat_soa1/carma_model_mod.F90 +++ b/src/physics/carma/models/trop_strat_soa1/carma_model_mod.F90 @@ -151,18 +151,6 @@ module carma_model_mod real(kind=f), public, parameter :: vmrat_MXAER = 2.2588_f !2.4610_f ! volume ratio ! Physics buffer index for sulfate surface area density - integer :: ipbuf4sadsulf = -1 ! total aerosol surface area density over all groups (cm2/cm3) - integer :: ipbuf4wtp = -1 ! sulfate weight percent H2SO4 composition - integer :: ipbuf4reffaer = -1 ! total aerosol effective radius over all groups (cm) - integer :: ipbuf4reff(NGROUP) = -1 ! aerosol effective radius per group (m) - integer :: ipbuf4numnkg(NGROUP) = -1 ! aerosol number density (#/kg air) - integer :: ipbuf4sad(NBIN,NGROUP) = -1 ! aerosol surface area density per bin (cm2/cm3) - integer :: ipbuf4elem1mr(NBIN,NGROUP) = -1 - integer :: ipbuf4binnkg(NBIN,NGROUP) = -1 - integer :: ipbuf4kappa(NBIN,NGROUP) = -1 ! hygroscopicity factor per bin - integer :: ipbuf4wetr(NBIN,NGROUP) = -1 ! aerosol wet radius per bin - integer :: ipbuf4dryr(NBIN,NGROUP) = -1 ! aerosol dry radius per bin - integer :: ipbuf4rmass(NBIN,NGROUP) = -1 ! aerosol mass per bin integer :: ipbuf4soa(NBIN) = -1 integer :: ipbuf4soacm(NBIN) = -1 integer :: ipbuf4soapt(NBIN) = -1 @@ -328,11 +316,11 @@ subroutine CARMAMODEL_DefineModel(carma, rc) ! should be 6 characters or less and without spaces. refidx(:,1) = CMPLX(shellreal(:), shellimag(:), kind=f) call CARMAELEMENT_Create(carma, I_ELEM_PRSUL, I_GRP_PRSUL, "Sulfate", & - RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, shortname="PRSUL", refidx=refidx) + RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, shortname="PRSULF", refidx=refidx) if (rc < 0) call endrun('CARMA_DefineModel::CARMA_AddElement failed.') call CARMAELEMENT_Create(carma, I_ELEM_MXAER, I_GRP_MXAER, "Sulfate in mixed sulfate", & - RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, kappa=Kappa_SULF, shortname="MXAER", refidx=refidx) + RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, kappa=Kappa_SULF, shortname="MXSULF", refidx=refidx) if (rc < 0) call endrun('CARMA_DefineModel::CARMA_AddElement failed.') call CARMAELEMENT_Create(carma, I_ELEM_MXOC, I_GRP_MXAER, "organic carbon", & @@ -400,31 +388,9 @@ subroutine CARMAMODEL_DefineModel(carma, rc) if (rc < 0) call endrun('carma_register::CARMAGROUP_Get failed.') !write(*,*) "igroup",igroup,"sname",sname - ! effective radius. ERMIXAER - call pbuf_add_field('ER'//trim(sname), 'global', dtype_r8, (/pcols, pver/), ipbuf4reff(igroup)) - !write(*,*) "'ER'//trim(sname)",'ER'//trim(sname) - - ! number density #/kg NBMXAER - call pbuf_add_field('NB'//trim(sname), 'global', dtype_r8, (/pcols, pver/), ipbuf4numnkg(igroup)) - !write(*,*) "'NB'//trim(sname)",'NB'//trim(sname) - ! sulfate mass and number density for each bin ! e.g. CRSULF01 first element mass mixing ratio; NBMXAER01 #/kg do ibin=1,NBIN - write (outputbin, "(I2.2)") ibin - write (outputname,"(A2)") trim(sname) - if (trim(sname)//outputbin .ne. outputname//"SULF"//outputbin) then - !write(*,*) "sname//outputbin",trim(sname)//outputbin - call pbuf_add_field(outputname//"SULF"//outputbin,'global', dtype_r8, (/pcols, pver/), ipbuf4elem1mr(ibin,igroup)) - !write(*,*) outputname//"SULF"//outputbin," ", trim(sname)//outputbin - end if - call pbuf_add_field("NB"//trim(sname)//outputbin,'global', dtype_r8, (/pcols, pver/), ipbuf4binnkg(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_kappa",'global', dtype_r8, (/pcols, pver/), ipbuf4kappa(ibin,igroup)) - !write(*,*) trim(sname)//outputbin//"_kappa" - call pbuf_add_field(trim(sname)//outputbin//"_wetr",'global', dtype_r8, (/pcols, pver/), ipbuf4wetr(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_dryr",'global', dtype_r8, (/pcols, pver/), ipbuf4dryr(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_rmass",'global', dtype_r8, (/pcols/), ipbuf4rmass(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_sad", 'global', dtype_r8, (/pcols, pver/), ipbuf4sad(ibin,igroup)) if (igroup==I_GRP_MXAER) then call pbuf_add_field("DQDT_MXSOA"//outputbin,'global',dtype_r8,(/pcols,pver/), ipbuf4soa(ibin)) call pbuf_add_field("MXSOA"//outputbin//"CM",'physpkg',dtype_r8,(/pcols,pver/), ipbuf4soacm(ibin)) @@ -433,12 +399,6 @@ subroutine CARMAMODEL_DefineModel(carma, rc) end do end do - ! total surface area density (cm2/cm3) - call pbuf_add_field('SADSULF', 'global', dtype_r8, (/pcols, pver/), ipbuf4sadsulf) - ! total effective radius (cm) for REFF_AERO history output - call pbuf_add_field('REFFAER', 'global', dtype_r8, (/pcols, pver/), ipbuf4reffaer) - ! weight percent H2SO4 - call pbuf_add_field('WTP','global', dtype_r8, (/pcols, pver/), ipbuf4wtp) ! no2 photolysis rate constant (/sec) call pbuf_add_field('JNO2', 'global', dtype_r8, (/pcols,pver/), ipbuf4jno2) @@ -502,7 +462,6 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli real(r8), pointer, dimension(:,:) :: jno2_rate !! jno2 tendency due to gas-aerosol exchange kg/kg/s real(r8), pointer, dimension(:,:) :: soacm !! aerosol tendency due to gas-aerosol exchange kg/kg/s real(r8), pointer, dimension(:,:) :: soapt !! aerosol tendency due to no2 photolysis kg/kg/s - real(r8) :: mmr_total(cstate%f_NZ)!! mass mixing ratio of a group (kg/kg) real(r8) :: mmr_core(cstate%f_NZ)!! mass mixing ratio of the core (kg/kg) real(r8) :: mmr_soa(cstate%f_NZ) !! mass mixing ratio of soa element (kg/kg) real(r8) :: mmr(cstate%f_NZ) !! mass mixing ratio per bin (kg/kg) @@ -523,15 +482,9 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli do ibin = 1, NBIN - ! Get the mass of the concentration element. - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc) - if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMASTATE_GetBin failed.') - ! Iterate over the core elements, looking for the SOA element. Once found, ! determine the new SOA taking into account both the addition of condensed ! SOA and the loss of photolyzed SOA. - ! - ! Both the total mass and the SOA mass need to be adjusted for the SOA change. do ielem = 1, ncore call CARMASTATE_GetBin(cstate, icorelem(ielem), ibin, mmr(:), rc) @@ -563,16 +516,10 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new values for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') - call CARMASTATE_SetBin(cstate, ienconc, ibin, mmr_total, rc) - if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') - exit end if !mxsoa end do !ielem @@ -614,7 +561,6 @@ subroutine CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state, pbuf, ptend, i real(r8) :: totreff(cstate%f_NZ) !! total volume density, used to calculate total effective radius (cm) for history output real(r8) :: reff(cstate%f_NZ) !! wet effective radius (m) real(r8) :: mmr(cstate%f_NZ) !! mass mixing ratio per bin (kg/kg) - real(r8) :: mmr_total(cstate%f_NZ)!! mass mixing ratio of a group (kg/kg) real(r8) :: coremmr(cstate%f_NZ) !! mmr of all the core real(r8) :: mmr_gas(cstate%f_NZ) !! gas mass mixing ratio (kg/kg) real(r8) :: numnkg(cstate%f_NZ) !! total number density (#/kg) @@ -633,153 +579,19 @@ subroutine CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state, pbuf, ptend, i real(r8), pointer, dimension(:,:) :: sadsulf_ptr !! Total surface area density pointer (cm2/cm3) real(r8), pointer, dimension(:,:) :: reffaer_ptr !! Total effective radius pointer (cm) for history output - real(r8), pointer, dimension(:,:) :: wtp_ptr !! weight percent pointer + real(r8), pointer, dimension(:,:) :: wtp_ptr !! weight percent pointer real(r8), pointer, dimension(:,:) :: sad_ptr !! Surface area density pointer real(r8), pointer, dimension(:,:) :: reff_ptr !! Effective radius pointer real(r8), pointer, dimension(:,:) :: numnkg_ptr !! Each group number density pointer real(r8), pointer, dimension(:,:) :: binnkg_ptr !! Each bin number density pointer real(r8), pointer, dimension(:,:) :: elem1mr_ptr !! First element mmr pointer - real(r8), pointer, dimension(:,:) :: kappa_ptr !! kappa pointer - real(r8), pointer, dimension(:,:) :: wetr_ptr !! wet radius pointer - real(r8), pointer, dimension(:,:) :: dryr_ptr !! dry radius + real(r8), pointer, dimension(:,:) :: kappa_ptr !! kappa pointer + real(r8), pointer, dimension(:,:) :: wetr_ptr !! wet radius pointer + real(r8), pointer, dimension(:,:) :: dryr_ptr !! dry radius ! Default return code. rc = RC_OK - ! Get the air density - call CARMASTATE_GetState(cstate, rc, rhoa_wet=rhoa_wet) - - totad(:) = 0.0_r8 ! total aerosol surface area density (cm2/cm3) - totreff(:) = 0.0_r8 ! total volume density, used to calculate total effective radius (cm) for REFF_AERO history output - - ! calculated SAD, RE, and number density (#/kg) for each group - do igroup = 1, NGROUP - - ad(:) = 0.0_r8 ! wet aerosol surface area density (cm2/cm3) - reff(:) = 0.0_r8 ! effective radius (m) - numnkg(:) = 0.0_r8 ! number density (#/kg) - - call CARMAGROUP_Get(carma, igroup, rc, ienconc=ienconc,ncore=ncore,icorelem=icorelem, rmass=rmass) - do ibin = 1, NBIN - - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, & - numberDensity=numberDensity, r_wet=r_wet) - if (rc < 0) call endrun('CARMA_DiagnoseBulk::CARMASTATE_GetBin failed.') - - if (ncore .ne. 0)then - do icore = 1,ncore - call CARMASTATE_GetBin(cstate, icorelem(icore), ibin, mmr(:), rc) - call CARMAELEMENT_Get(carma, icorelem(icore), rc, icomposition=icomposition) - end do - end if - - ! Calculate the total densities. - ! NOTE: Calculate AD in cm2/cm3. - if (numberDensity(1) /= CAM_FILL) then - ad(:) = ad(:) + numberDensity(:) * (r_wet(:)**2) - reff(:) = reff(:) + numberDensity(:) * (r_wet(:)**3) - end if - end do - - totreff(:) = totreff(:) + reff(:) ! volume density in cm3 - reff(:) = reff(:) / ad(:) ! wet effective radius in cm - reff(:) = reff(:) / 100.0_r8 ! cm -> m - - ad(:) = ad(:) * 4.0_r8 * PI ! surface area density in cm2/cm3 - - ! airdensity from carma state - ! convert the number density from #/cm3 to #/kg - ! number Density #/cm3; rhoa_wet kg/m3 - numnkg(:) = 1.e6_r8*numberDensity(:)/rhoa_wet(:) !#/kg - - call pbuf_get_field(pbuf, ipbuf4reff(igroup), reff_ptr) ! effective radius m - call pbuf_get_field(pbuf, ipbuf4numnkg(igroup), numnkg_ptr) ! number density #/kg - - ! put variables in physics buffer - reff_ptr(icol, :cstate%f_NZ) = reff(:cstate%f_NZ) - numnkg_ptr(icol, :cstate%f_NZ)= numnkg(:cstate%f_NZ) - - !calculate the surface area density including all bins (cm2/cm3) - totad(:) = totad(:)+ad(:) - - end do - - totreff(:) = 4.0_r8 * PI * totreff(:) / totad(:) ! cm - - call pbuf_get_field(pbuf, ipbuf4sadsulf, sadsulf_ptr) ! surface area density (cm2/cm3) - sadsulf_ptr(icol, :cstate%f_NZ) = totad(:cstate%f_NZ) - - call pbuf_get_field(pbuf, ipbuf4reffaer, reffaer_ptr) ! effective radius (cm) - reffaer_ptr(icol, :cstate%f_NZ) = totreff(:cstate%f_NZ) - - do igas = 1,NGAS - if(igas .eq. I_GAS_H2SO4)then ! only output the sulfate weight percent - call CARMASTATE_GetGas(cstate, igas, mmr_gas(:), rc, wtpct=wtpct) - call pbuf_get_field(pbuf, ipbuf4wtp, wtp_ptr) - wtp_ptr(icol, :cstate%f_NZ) = wtpct(:cstate%f_NZ) - end if - end do - - ! add new CRSULF element that is only the sulfur part from MIXAER: - ! calculate CRSULF mass mixing ratio : MXAER minus CROC+CRBC+CRBC+CRDUST+CRSALT (+SUM(CRSOA*) - - do igroup = 1, NGROUP - call CARMAGROUP_Get(carma, igroup, rc, shortname=sname,ienconc=ienconc,ncore=ncore,icorelem=icorelem, rmass=rmass) - !write(*,*) "igroup",igroup,"ncore",ncore,"ienconc",ienconc,"icorelem",icorelem - do ibin = 1, NBIN - - elem1mr(:) = 0._r8 - binnkg(:) = 0._r8 - coremmr(:) = 0._r8 - - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, numberDensity=numberDensity, kappa=kappa, r_wet=r_wet,rhop_dry=rhop_dry) - - if (ncore .ne. 0)then - do icore = 1,ncore - call CARMASTATE_GetBin(cstate, icorelem(icore), ibin, mmr(:), rc) - coremmr(:) = coremmr(:) + mmr(:) - end do - - if (any(coremmr(:) .gt. mmr_total(:))) then - where(coremmr(:) .gt. mmr_total(:)) mmr_total = coremmr - call CARMASTATE_SetBin(cstate, ienconc, ibin, mmr_total(:), rc) - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, numberDensity=numberDensity, r_wet=r_wet,rhop_dry=rhop_dry) - end if - elem1mr(:) = mmr_total(:)-coremmr(:) - elem1mr(:) = max(elem1mr(:),0._f) - - call pbuf_get_field(pbuf, ipbuf4elem1mr(ibin,igroup), elem1mr_ptr) - elem1mr_ptr(icol, :cstate%f_NZ) = elem1mr(:cstate%f_NZ) - else - elem1mr(:) = mmr_total(:) - call pbuf_get_field(pbuf, ipbuf4elem1mr(ibin,igroup), elem1mr_ptr) - elem1mr_ptr(icol, :cstate%f_NZ) = elem1mr(:cstate%f_NZ) - end if - binnkg(:) = 1.e6_r8*numberDensity(:)/rhoa_wet(:) !#/kg - - call pbuf_get_field(pbuf, ipbuf4binnkg(ibin,igroup), binnkg_ptr) - call pbuf_get_field(pbuf, ipbuf4kappa(ibin,igroup), kappa_ptr) - - binnkg_ptr(icol, :cstate%f_NZ) = binnkg(:cstate%f_NZ) - kappa_ptr(icol, :cstate%f_NZ) = kappa(:cstate%f_NZ) - - call pbuf_get_field(pbuf, ipbuf4wetr(ibin,igroup), wetr_ptr) - call pbuf_get_field(pbuf, ipbuf4dryr(ibin,igroup), dryr_ptr) - call pbuf_get_field(pbuf, ipbuf4sad(ibin,igroup), sad_ptr) - - wetr_ptr(icol, :cstate%f_NZ) = r_wet(:cstate%f_NZ) - !dryr_ptr(icol, :cstate%f_NZ) = (rmass(ibin)/(4._f/3._f*PI*rhop_dry(:cstate%f_NZ)))**(1._f/3._f) !cm - ! r = (mass / rho / 4 * 3 / PI)^(1/3) - dryr_ptr(icol, :cstate%f_NZ) = (rmass(ibin)/rhop_dry(:cstate%f_NZ) / 4._f * 3._f / PI)**(1._f/3._f) !cm - - sad_ptr(icol, :cstate%f_NZ) = 4.0_r8 * PI * numberDensity(:cstate%f_NZ) * (r_wet(:cstate%f_NZ)**2) !cm2/cm3 - - !write(*,*) 'CARMA igroup, ibin, rmass(ibin), rhop_dry',igroup, ibin, rmass(ibin), rhop_dry(:cstate%f_NZ) - !write(*,*) 'CARMA dryr igroup, ibin',igroup, ibin, dryr_ptr(icol, :cstate%f_NZ) - - end do !NBIN - end do !NGROUP - return end subroutine CARMAMODEL_DiagnoseBulk @@ -926,7 +738,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Organic carbon emssions - if ((ielem == I_ELEM_MXOC) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXOC) then if (carma_BCOCemissions == 'Yu2015') then call get_lat_all_p(lchnk, ncol, ilat) call get_lon_all_p(lchnk, ncol, ilon) @@ -942,7 +754,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Black carbon emissions - if ((ielem == I_ELEM_MXBC) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXBC) then if (carma_BCOCemissions == 'Yu2015') then do icol = 1,ncol smoke(icol) = BCnew(ilat(icol), ilon(icol), mon) @@ -991,7 +803,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Dust emissions - if ((ielem == I_ELEM_MXDUST) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXDUST) then ! The radius should be determined by the dust density not the group ! density @@ -1044,7 +856,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, ! Sea salt emissions - if ((ielem == I_ELEM_MXSALT) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXSALT) then ! The radius should be determined by the dust density not the group ! density @@ -1216,33 +1028,16 @@ subroutine CARMAMODEL_InitializeModel(carma, lq_carma, pbuf2d, rc) end if if (is_first_step()) then + ! Initialize physics buffer fields do igroup = 1, NGROUP - call pbuf_set_field(pbuf2d, ipbuf4reff(igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4numnkg(igroup), 0.0_r8 ) - - call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) - if (RC /= RC_OK) return - - do ibin=1,NBIN - if (ipbuf4elem1mr(ibin,igroup)>0) then - call pbuf_set_field(pbuf2d, ipbuf4elem1mr(ibin,igroup), 0.0_r8 ) - endif - call pbuf_set_field(pbuf2d, ipbuf4binnkg(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4kappa(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4wetr(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4dryr(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4sad(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4rmass(ibin,igroup), 1.0e-3_r8*rmass(ibin) ) ! convert rmass from g to kg + do ibin = 1, NBIN if (igroup==I_GRP_MXAER) then call pbuf_set_field(pbuf2d, ipbuf4soa(ibin), 0.0_r8 ) end if end do end do - call pbuf_set_field(pbuf2d, ipbuf4sadsulf, 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4reffaer, 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4wtp, 0.0_r8 ) call pbuf_set_field(pbuf2d, ipbuf4jno2, 0.0_r8 ) endif @@ -2896,8 +2691,6 @@ subroutine CARMAMODEL_CalculateCloudborneDiagnostics(carma, state, pbuf, aercldd else if (shortname .eq. "MXSOA") then bdsoa(:ncols, :) = bdsoa(:ncols, :) + mmr(:ncols,:) * mair(:ncols,:) end if - - mixso4(:ncols,:) = mixso4(:ncols,:) - mmr(:ncols,:) * mair(:ncols,:) end do end do @@ -3011,8 +2804,6 @@ subroutine CARMAMODEL_OutputBudgetDiagnostics(carma, icnst4elem, icnst4gas, stat if (ptend%lq(icnst)) then tottend(:) = tottend(:) - ptend%q(icol,:,icnst) * mair(:) end if - bdmxso4(icol) = bdmxso4(icol) - sum(state%q(icol,:,icnst) * mair(:)) - cmxflux(icol) = cmxflux(icol) - (cflux(icol,icnst) - old_cflux(icol,icnst)) end do mixtend(icol) = mixtend(icol) + sum(tottend(:)) @@ -3225,9 +3016,6 @@ subroutine CARMAMODEL_OutputDiagnostics(carma, icnst4elem, state, ptend, pbuf, c do i = 1, ncore icnst = icnst4elem(icorelem(i), ibin) - mixso4mr(icol,:) = mixso4mr(icol,:) - state%q(icol,:,icnst) - mixso4(icol) = mixso4(icol) - sum(state%q(icol,:,icnst) * mair(:)) - cmxflux(icol) = cmxflux(icol) - cam_in%cflx(icol,icnst) call CARMAELEMENT_Get(carma, icorelem(i), rc, shortname=shortname) if (shortname .eq. "MXBC") then @@ -3928,7 +3716,7 @@ end subroutine CARMAMODEL_SurfaceWind !! could use /components/cam/src/chemistry/aerosol/soil_erod_mod.F90 here insted of this routine? subroutine CARMAMODEL_ReadSoilErosionFactor(rc) use ppgrid, only: begchunk, endchunk, pcols - use ioFileMod, only: getfil + use ioFileMod, only: getfil use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish use phys_grid, only: get_rlon_all_p, get_rlat_all_p, get_ncols_p use wrap_nf diff --git a/src/physics/carma/models/trop_strat_soa5/carma_model_mod.F90 b/src/physics/carma/models/trop_strat_soa5/carma_model_mod.F90 index 43e29455dc..4b74d7f2c5 100644 --- a/src/physics/carma/models/trop_strat_soa5/carma_model_mod.F90 +++ b/src/physics/carma/models/trop_strat_soa5/carma_model_mod.F90 @@ -159,18 +159,6 @@ module carma_model_mod real(kind=f), public, parameter :: vmrat_MXAER = 2.2588_f !2.4610_f ! volume ratio ! Physics buffer index for sulfate surface area density - integer :: ipbuf4sadsulf = -1 ! total aerosol surface area density over all groups (cm2/cm3) - integer :: ipbuf4wtp = -1 ! sulfate weight percent H2SO4 composition - integer :: ipbuf4reffaer = -1 ! total aerosol effective radius over all groups (cm) - integer :: ipbuf4reff(NGROUP) = -1 ! aerosol effective radius per group (m) - integer :: ipbuf4numnkg(NGROUP) = -1 ! aerosol number density (#/kg air) - integer :: ipbuf4sad(NBIN,NGROUP) = -1 ! aerosol surface area density per bin (cm2/cm3) - integer :: ipbuf4elem1mr(NBIN,NGROUP) = -1 - integer :: ipbuf4binnkg(NBIN,NGROUP) = -1 - integer :: ipbuf4kappa(NBIN,NGROUP) = -1 ! hygroscopicity factor per bin - integer :: ipbuf4wetr(NBIN,NGROUP) = -1 ! aerosol wet radius per bin - integer :: ipbuf4dryr(NBIN,NGROUP) = -1 ! aerosol dry radius per bin - integer :: ipbuf4rmass(NBIN,NGROUP) = -1 ! aerosol mass per bin integer :: ipbuf4soa1(NBIN) = -1 integer :: ipbuf4soa2(NBIN) = -1 integer :: ipbuf4soa3(NBIN) = -1 @@ -348,11 +336,11 @@ subroutine CARMAMODEL_DefineModel(carma, rc) ! should be 6 characters or less and without spaces. refidx(:,1) = CMPLX(shellreal(:), shellimag(:), kind=f) call CARMAELEMENT_Create(carma, I_ELEM_PRSUL, I_GRP_PRSUL, "Sulfate", & - RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, shortname="PRSUL", refidx=refidx) + RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, shortname="PRSULF", refidx=refidx) if (rc < 0) call endrun('CARMA_DefineModel::CARMA_AddElement failed.') call CARMAELEMENT_Create(carma, I_ELEM_MXAER, I_GRP_MXAER, "Sulfate in mixed sulfate", & - RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, kappa=Kappa_SULF, shortname="MXAER", refidx=refidx) + RHO_SULFATE, I_VOLATILE, I_H2SO4, rc, kappa=Kappa_SULF, shortname="MXSULF", refidx=refidx) if (rc < 0) call endrun('CARMA_DefineModel::CARMA_AddElement failed.') call CARMAELEMENT_Create(carma, I_ELEM_MXOC, I_GRP_MXAER, "organic carbon", & @@ -436,31 +424,9 @@ subroutine CARMAMODEL_DefineModel(carma, rc) if (rc < 0) call endrun('carma_register::CARMAGROUP_Get failed.') !write(*,*) "igroup",igroup,"sname",sname - ! effective radius. ERMIXAER - call pbuf_add_field('ER'//trim(sname), 'global', dtype_r8, (/pcols, pver/), ipbuf4reff(igroup)) - !write(*,*) "'ER'//trim(sname)",'ER'//trim(sname) - - ! number density #/kg NBMXAER - call pbuf_add_field('NB'//trim(sname), 'global', dtype_r8, (/pcols, pver/), ipbuf4numnkg(igroup)) - !write(*,*) "'NB'//trim(sname)",'NB'//trim(sname) - ! sulfate mass and number density for each bin ! e.g. CRSULF01 first element mass mixing ratio; NBMXAER01 #/kg do ibin=1,NBIN - write (outputbin, "(I2.2)") ibin - write (outputname,"(A2)") trim(sname) - if (trim(sname)//outputbin .ne. outputname//"SULF"//outputbin) then - !write(*,*) "sname//outputbin",trim(sname)//outputbin - call pbuf_add_field(outputname//"SULF"//outputbin,'global', dtype_r8, (/pcols, pver/), ipbuf4elem1mr(ibin,igroup)) - !write(*,*) outputname//"SULF"//outputbin," ", trim(sname)//outputbin - end if - call pbuf_add_field("NB"//trim(sname)//outputbin,'global', dtype_r8, (/pcols, pver/), ipbuf4binnkg(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_kappa",'global', dtype_r8, (/pcols, pver/), ipbuf4kappa(ibin,igroup)) - !write(*,*) trim(sname)//outputbin//"_kappa" - call pbuf_add_field(trim(sname)//outputbin//"_wetr",'global', dtype_r8, (/pcols, pver/), ipbuf4wetr(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_dryr",'global', dtype_r8, (/pcols, pver/), ipbuf4dryr(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_rmass",'global', dtype_r8, (/pcols/), ipbuf4rmass(ibin,igroup)) - call pbuf_add_field(trim(sname)//outputbin//"_sad", 'global', dtype_r8, (/pcols, pver/), ipbuf4sad(ibin,igroup)) if (igroup==I_GRP_MXAER) then call pbuf_add_field("DQDT_MXSOA1"//outputbin,'global',dtype_r8,(/pcols,pver/), ipbuf4soa1(ibin)) call pbuf_add_field("DQDT_MXSOA2"//outputbin,'global',dtype_r8,(/pcols,pver/), ipbuf4soa2(ibin)) @@ -481,12 +447,6 @@ subroutine CARMAMODEL_DefineModel(carma, rc) end do end do - ! total surface area density (cm2/cm3) - call pbuf_add_field('SADSULF', 'global', dtype_r8, (/pcols, pver/), ipbuf4sadsulf) - ! total effective radius (cm) for REFF_AERO history output - call pbuf_add_field('REFFAER', 'global', dtype_r8, (/pcols, pver/), ipbuf4reffaer) - ! weight percent H2SO4 - call pbuf_add_field('WTP','global', dtype_r8, (/pcols, pver/), ipbuf4wtp) ! no2 photolysis rate constant (/sec) call pbuf_add_field('JNO2', 'global', dtype_r8, (/pcols,pver/), ipbuf4jno2) @@ -558,7 +518,6 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli real(r8), pointer, dimension(:,:) :: soapt3 !! aerosol tendency due to no2 photolysis kg/kg/s real(r8), pointer, dimension(:,:) :: soapt4 !! aerosol tendency due to no2 photolysis kg/kg/s real(r8), pointer, dimension(:,:) :: soapt5 !! aerosol tendency due to no2 photolysis kg/kg/s - real(r8) :: mmr_total(cstate%f_NZ)!! mass mixing ratio of a group (kg/kg) real(r8) :: mmr_core(cstate%f_NZ)!! mass mixing ratio of the core (kg/kg) real(r8) :: mmr_soa(cstate%f_NZ) !! mass mixing ratio of soa element (kg/kg) real(r8) :: mmr(cstate%f_NZ) !! mass mixing ratio per bin (kg/kg) @@ -579,15 +538,9 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli do ibin = 1, NBIN - ! Get the mass of the concentration element. - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc) - if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMASTATE_GetBin failed.') - - ! Iterate over the core elements, looking for the SOA element. Once found, + ! Iterate over the core elements, looking for the SOA elements. Once found, ! determine the new SOA taking into account both the addition of condensed ! SOA and the loss of photolyzed SOA. - ! - ! Both the total mass and the SOA mass need to be adjusted for the SOA change. do ielem = 1, ncore call CARMASTATE_GetBin(cstate, icorelem(ielem), ibin, mmr(:), rc) @@ -619,10 +572,7 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt1(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new value for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') @@ -650,10 +600,7 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt2(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new value for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') @@ -681,10 +628,7 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt3(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new value for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') @@ -712,10 +656,7 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt4(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new value for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') @@ -743,19 +684,12 @@ subroutine CARMAMODEL_DiagnoseBins(carma, cstate, state, pbuf, icol, dt, rc, rli ! chemistry. mmr_soa(:) = max(0.0_r8, mmr_soa(:) + soapt5(icol,:) * dt) - ! Adjust the total MMR. - mmr_total(:) = mmr_total(:) + (mmr_soa(:) - mmr(:)) - - ! Save out these new values for SOA and total MMR. + ! Save out these new value for SOA. call CARMASTATE_SetBin(cstate, icorelem(ielem), ibin, mmr_soa, rc) if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') end if !mxsoa5 end do !ielem - - call CARMASTATE_SetBin(cstate, ienconc, ibin, mmr_total, rc) - if (rc /= RC_OK) call endrun('CARMA_DiagnoseBins::CARMAGROUP_SetBin failed.') - end do !ibin end subroutine CARMAMODEL_DiagnoseBins @@ -794,7 +728,6 @@ subroutine CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state, pbuf, ptend, i real(r8) :: totreff(cstate%f_NZ) !! total volume density, used to calculate total effective radius (cm) for history output real(r8) :: reff(cstate%f_NZ) !! wet effective radius (m) real(r8) :: mmr(cstate%f_NZ) !! mass mixing ratio per bin (kg/kg) - real(r8) :: mmr_total(cstate%f_NZ)!! mass mixing ratio of a group (kg/kg) real(r8) :: coremmr(cstate%f_NZ) !! mmr of all the core real(r8) :: mmr_gas(cstate%f_NZ) !! gas mass mixing ratio (kg/kg) real(r8) :: numnkg(cstate%f_NZ) !! total number density (#/kg) @@ -813,153 +746,19 @@ subroutine CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state, pbuf, ptend, i real(r8), pointer, dimension(:,:) :: sadsulf_ptr !! Total surface area density pointer (cm2/cm3) real(r8), pointer, dimension(:,:) :: reffaer_ptr !! Total effective radius pointer (cm) for history output - real(r8), pointer, dimension(:,:) :: wtp_ptr !! weight percent pointer + real(r8), pointer, dimension(:,:) :: wtp_ptr !! weight percent pointer real(r8), pointer, dimension(:,:) :: sad_ptr !! Surface area density pointer real(r8), pointer, dimension(:,:) :: reff_ptr !! Effective radius pointer real(r8), pointer, dimension(:,:) :: numnkg_ptr !! Each group number density pointer real(r8), pointer, dimension(:,:) :: binnkg_ptr !! Each bin number density pointer real(r8), pointer, dimension(:,:) :: elem1mr_ptr !! First element mmr pointer - real(r8), pointer, dimension(:,:) :: kappa_ptr !! kappa pointer - real(r8), pointer, dimension(:,:) :: wetr_ptr !! wet radius pointer - real(r8), pointer, dimension(:,:) :: dryr_ptr !! dry radius + real(r8), pointer, dimension(:,:) :: kappa_ptr !! kappa pointer + real(r8), pointer, dimension(:,:) :: wetr_ptr !! wet radius pointer + real(r8), pointer, dimension(:,:) :: dryr_ptr !! dry radius ! Default return code. rc = RC_OK - ! Get the air density - call CARMASTATE_GetState(cstate, rc, rhoa_wet=rhoa_wet) - - totad(:) = 0.0_r8 ! total aerosol surface area density (cm2/cm3) - totreff(:) = 0.0_r8 ! total volume density, used to calculate total effective radius (cm) for REFF_AERO history output - - ! calculated SAD, RE, and number density (#/kg) for each group - do igroup = 1, NGROUP - - ad(:) = 0.0_r8 ! wet aerosol surface area density (cm2/cm3) - reff(:) = 0.0_r8 ! effective radius (m) - numnkg(:) = 0.0_r8 ! number density (#/kg) - - call CARMAGROUP_Get(carma, igroup, rc, ienconc=ienconc,ncore=ncore,icorelem=icorelem, rmass=rmass) - do ibin = 1, NBIN - - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, & - numberDensity=numberDensity, r_wet=r_wet) - if (rc < 0) call endrun('CARMA_DiagnoseBulk::CARMASTATE_GetBin failed.') - - if (ncore .ne. 0)then - do icore = 1,ncore - call CARMASTATE_GetBin(cstate, icorelem(icore), ibin, mmr(:), rc) - call CARMAELEMENT_Get(carma, icorelem(icore), rc, icomposition=icomposition) - end do - end if - - ! Calculate the total densities. - ! NOTE: Calculate AD in cm2/cm3. - if (numberDensity(1) /= CAM_FILL) then - ad(:) = ad(:) + numberDensity(:) * (r_wet(:)**2) - reff(:) = reff(:) + numberDensity(:) * (r_wet(:)**3) - end if - end do - - totreff(:) = totreff(:) + reff(:) ! volume density in cm3 - reff(:) = reff(:) / ad(:) ! wet effective radius in cm - reff(:) = reff(:) / 100.0_r8 ! cm -> m - - ad(:) = ad(:) * 4.0_r8 * PI ! surface area density in cm2/cm3 - - ! airdensity from carma state - ! convert the number density from #/cm3 to #/kg - ! number Density #/cm3; rhoa_wet kg/m3 - numnkg(:) = 1.e6_r8*numberDensity(:)/rhoa_wet(:) !#/kg - - call pbuf_get_field(pbuf, ipbuf4reff(igroup), reff_ptr) ! effective radius m - call pbuf_get_field(pbuf, ipbuf4numnkg(igroup), numnkg_ptr) ! number density #/kg - - ! put variables in physics buffer - reff_ptr(icol, :cstate%f_NZ) = reff(:cstate%f_NZ) - numnkg_ptr(icol, :cstate%f_NZ)= numnkg(:cstate%f_NZ) - - !calculate the surface area density including all bins (cm2/cm3) - totad(:) = totad(:)+ad(:) - - end do - - totreff(:) = 4.0_r8 * PI * totreff(:) / totad(:) ! cm - - call pbuf_get_field(pbuf, ipbuf4sadsulf, sadsulf_ptr) ! surface area density (cm2/cm3) - sadsulf_ptr(icol, :cstate%f_NZ) = totad(:cstate%f_NZ) - - call pbuf_get_field(pbuf, ipbuf4reffaer, reffaer_ptr) ! effective radius (cm) - reffaer_ptr(icol, :cstate%f_NZ) = totreff(:cstate%f_NZ) - - do igas = 1,NGAS - if(igas .eq. I_GAS_H2SO4)then ! only output the sulfate weight percent - call CARMASTATE_GetGas(cstate, igas, mmr_gas(:), rc, wtpct=wtpct) - call pbuf_get_field(pbuf, ipbuf4wtp, wtp_ptr) - wtp_ptr(icol, :cstate%f_NZ) = wtpct(:cstate%f_NZ) - end if - end do - - ! add new CRSULF element that is only the sulfur part from MIXAER: - ! calculate CRSULF mass mixing ratio : MXAER minus CROC+CRBC+CRBC+CRDUST+CRSALT (+SUM(CRSOA*) - - do igroup = 1, NGROUP - call CARMAGROUP_Get(carma, igroup, rc, shortname=sname,ienconc=ienconc,ncore=ncore,icorelem=icorelem, rmass=rmass) - !write(*,*) "igroup",igroup,"ncore",ncore,"ienconc",ienconc,"icorelem",icorelem - do ibin = 1, NBIN - - elem1mr(:) = 0._r8 - binnkg(:) = 0._r8 - coremmr(:) = 0._r8 - - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, numberDensity=numberDensity, kappa=kappa, r_wet=r_wet,rhop_dry=rhop_dry) - - if (ncore .ne. 0)then - do icore = 1,ncore - call CARMASTATE_GetBin(cstate, icorelem(icore), ibin, mmr(:), rc) - coremmr(:) = coremmr(:) + mmr(:) - end do - - if (any(coremmr(:) .gt. mmr_total(:))) then - where(coremmr(:) .gt. mmr_total(:)) mmr_total = coremmr - call CARMASTATE_SetBin(cstate, ienconc, ibin, mmr_total(:), rc) - call CARMASTATE_GetBin(cstate, ienconc, ibin, mmr_total(:), rc, numberDensity=numberDensity, r_wet=r_wet,rhop_dry=rhop_dry) - end if - elem1mr(:) = mmr_total(:)-coremmr(:) - elem1mr(:) = max(elem1mr(:),0._f) - - call pbuf_get_field(pbuf, ipbuf4elem1mr(ibin,igroup), elem1mr_ptr) - elem1mr_ptr(icol, :cstate%f_NZ) = elem1mr(:cstate%f_NZ) - else - elem1mr(:) = mmr_total(:) - call pbuf_get_field(pbuf, ipbuf4elem1mr(ibin,igroup), elem1mr_ptr) - elem1mr_ptr(icol, :cstate%f_NZ) = elem1mr(:cstate%f_NZ) - end if - binnkg(:) = 1.e6_r8*numberDensity(:)/rhoa_wet(:) !#/kg - - call pbuf_get_field(pbuf, ipbuf4binnkg(ibin,igroup), binnkg_ptr) - call pbuf_get_field(pbuf, ipbuf4kappa(ibin,igroup), kappa_ptr) - - binnkg_ptr(icol, :cstate%f_NZ) = binnkg(:cstate%f_NZ) - kappa_ptr(icol, :cstate%f_NZ) = kappa(:cstate%f_NZ) - - call pbuf_get_field(pbuf, ipbuf4wetr(ibin,igroup), wetr_ptr) - call pbuf_get_field(pbuf, ipbuf4dryr(ibin,igroup), dryr_ptr) - call pbuf_get_field(pbuf, ipbuf4sad(ibin,igroup), sad_ptr) - - wetr_ptr(icol, :cstate%f_NZ) = r_wet(:cstate%f_NZ) - !dryr_ptr(icol, :cstate%f_NZ) = (rmass(ibin)/(4._f/3._f*PI*rhop_dry(:cstate%f_NZ)))**(1._f/3._f) !cm - ! r = (mass / rho / 4 * 3 / PI)^(1/3) - dryr_ptr(icol, :cstate%f_NZ) = (rmass(ibin)/rhop_dry(:cstate%f_NZ) / 4._f * 3._f / PI)**(1._f/3._f) !cm - - sad_ptr(icol, :cstate%f_NZ) = 4.0_r8 * PI * numberDensity(:cstate%f_NZ) * (r_wet(:cstate%f_NZ)**2) !cm2/cm3 - - !write(*,*) 'CARMA igroup, ibin, rmass(ibin), rhop_dry',igroup, ibin, rmass(ibin), rhop_dry(:cstate%f_NZ) - !write(*,*) 'CARMA dryr igroup, ibin',igroup, ibin, dryr_ptr(icol, :cstate%f_NZ) - - end do !NBIN - end do !NGROUP - return end subroutine CARMAMODEL_DiagnoseBulk @@ -1106,7 +905,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Organic carbon emssions - if ((ielem == I_ELEM_MXOC) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXOC) then if (carma_BCOCemissions == 'Yu2015') then call get_lat_all_p(lchnk, ncol, ilat) call get_lon_all_p(lchnk, ncol, ilon) @@ -1122,7 +921,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Black carbon emissions - if ((ielem == I_ELEM_MXBC) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXBC) then if (carma_BCOCemissions == 'Yu2015') then do icol = 1,ncol smoke(icol) = BCnew(ilat(icol), ilon(icol), mon) @@ -1171,7 +970,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, end if ! Dust emissions - if ((ielem == I_ELEM_MXDUST) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXDUST) then ! The radius should be determined by the dust density not the group ! density @@ -1224,7 +1023,7 @@ subroutine CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, ! Sea salt emissions - if ((ielem == I_ELEM_MXSALT) .or. (ielem == I_ELEM_MXAER)) then + if (ielem == I_ELEM_MXSALT) then ! The radius should be determined by the dust density not the group ! density @@ -1396,24 +1195,10 @@ subroutine CARMAMODEL_InitializeModel(carma, lq_carma, pbuf2d, rc) end if if (is_first_step()) then + ! Initialize physics buffer fields do igroup = 1, NGROUP - call pbuf_set_field(pbuf2d, ipbuf4reff(igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4numnkg(igroup), 0.0_r8 ) - - call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) - if (RC /= RC_OK) return - - do ibin=1,NBIN - if (ipbuf4elem1mr(ibin,igroup)>0) then - call pbuf_set_field(pbuf2d, ipbuf4elem1mr(ibin,igroup), 0.0_r8 ) - endif - call pbuf_set_field(pbuf2d, ipbuf4binnkg(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4kappa(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4wetr(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4dryr(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4sad(ibin,igroup), 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4rmass(ibin,igroup), 1.0e-3_r8*rmass(ibin) ) ! convert rmass from g to kg + do ibin = 1, NBIN if (igroup==I_GRP_MXAER) then call pbuf_set_field(pbuf2d, ipbuf4soa1(ibin), 0.0_r8 ) call pbuf_set_field(pbuf2d, ipbuf4soa2(ibin), 0.0_r8 ) @@ -1424,9 +1209,6 @@ subroutine CARMAMODEL_InitializeModel(carma, lq_carma, pbuf2d, rc) end do end do - call pbuf_set_field(pbuf2d, ipbuf4sadsulf, 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4reffaer, 0.0_r8 ) - call pbuf_set_field(pbuf2d, ipbuf4wtp, 0.0_r8 ) call pbuf_set_field(pbuf2d, ipbuf4jno2, 0.0_r8 ) endif @@ -3123,8 +2905,6 @@ subroutine CARMAMODEL_CalculateCloudborneDiagnostics(carma, state, pbuf, aercldd else if (shortname .eq. "MXSOA5") then bdsoa5(:ncols, :) = bdsoa5(:ncols, :) + mmr(:ncols,:) * mair(:ncols,:) end if - - mixso4(:ncols,:) = mixso4(:ncols,:) - mmr(:ncols,:) * mair(:ncols,:) end do end do @@ -3242,8 +3022,6 @@ subroutine CARMAMODEL_OutputBudgetDiagnostics(carma, icnst4elem, icnst4gas, stat if (ptend%lq(icnst)) then tottend(:) = tottend(:) - ptend%q(icol,:,icnst) * mair(:) end if - bdmxso4(icol) = bdmxso4(icol) - sum(state%q(icol,:,icnst) * mair(:)) - cmxflux(icol) = cmxflux(icol) - (cflux(icol,icnst) - old_cflux(icol,icnst)) end do mixtend(icol) = mixtend(icol) + sum(tottend(:)) @@ -3499,9 +3277,6 @@ subroutine CARMAMODEL_OutputDiagnostics(carma, icnst4elem, state, ptend, pbuf, c do i = 1, ncore icnst = icnst4elem(icorelem(i), ibin) - mixso4mr(icol,:) = mixso4mr(icol,:) - state%q(icol,:,icnst) - mixso4(icol) = mixso4(icol) - sum(state%q(icol,:,icnst) * mair(:)) - cmxflux(icol) = cmxflux(icol) - cam_in%cflx(icol,icnst) call CARMAELEMENT_Get(carma, icorelem(i), rc, shortname=shortname) if (shortname .eq. "MXBC") then @@ -4213,7 +3988,7 @@ end subroutine CARMAMODEL_SurfaceWind !! could use /components/cam/src/chemistry/aerosol/soil_erod_mod.F90 here insted of this routine? subroutine CARMAMODEL_ReadSoilErosionFactor(rc) use ppgrid, only: begchunk, endchunk, pcols - use ioFileMod, only: getfil + use ioFileMod, only: getfil use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish use phys_grid, only: get_rlon_all_p, get_rlat_all_p, get_ncols_p use wrap_nf