From 44ed90f513dc870f8fe2d2286937417ddef4d3ad Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 14 Nov 2023 08:25:29 -0700 Subject: [PATCH] functional wghtpct routine modified: src/chemistry/aerosol/carma_aerosol_properties_mod.F90 modified: src/chemistry/aerosol/carma_aerosol_state_mod.F90 modified: src/physics/carma/cam/carma_intr.F90 --- .../aerosol/carma_aerosol_properties_mod.F90 | 22 ++-- .../aerosol/carma_aerosol_state_mod.F90 | 4 +- src/physics/carma/cam/carma_intr.F90 | 120 +++++++++++++----- 3 files changed, 105 insertions(+), 41 deletions(-) diff --git a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 index 4a0f45e9c8..ab56a49453 100644 --- a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 @@ -561,7 +561,8 @@ end function bin_name function scav_radius(self, bin_ndx) result(radius) use carma_model_mod, only: NBIN - use carma_intr, only :carma_get_bin_rmass + use carma_intr, only: carma_get_bin_rmass + use carma_intr, only: carma_get_group_by_name class(carma_aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin number @@ -570,15 +571,20 @@ function scav_radius(self, bin_ndx) result(radius) real(r8) :: mass ! the bin mass (g) real(r8) :: rho ! density (kg/m3) - integer :: igroup, ibin, rc, ispec + integer :: ispec character(len=32) :: spectype - ibin = bin_ndx - igroup = 1 - if (bin_ndx>NBIN) then - igroup = 2 - ibin = ibin-NBIN - end if + character(len=aero_name_len) :: bin_name, shortname + integer :: igroup, ibin, rc, nchr + + call rad_cnst_get_info_by_bin(0, bin_ndx, bin_name=bin_name) + + nchr = len_trim(bin_name)-2 + shortname = bin_name(:nchr) + + call carma_get_group_by_name(shortname, igroup, rc) + + read(bin_name(nchr+1:),*) ibin call carma_get_bin_rmass(igroup, ibin, mass, rc) diff --git a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 index b87d6bc68a..4f31e4cc6e 100644 --- a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 @@ -11,7 +11,7 @@ module carma_aerosol_state_mod use physconst, only: pi use carma_intr, only: carma_get_total_mmr, carma_get_dry_radius, carma_get_number, carma_get_number_cld use carma_intr, only: carma_get_group_by_name, carma_get_kappa, carma_get_dry_radius, carma_get_wet_radius - use carma_intr, only: carma_sulfate_wght_pct + use carma_intr, only: carma_get_wght_pct use ppgrid, only: begchunk, endchunk, pcols, pver implicit none @@ -452,7 +452,7 @@ function wgtpct(self, icol,ilev) result(wtp) integer, intent(in) :: icol,ilev real(r8) :: wtp ! weight precent of H2SO4/H2O solution for given icol, ilev - wtp = carma_sulfate_wght_pct(icol,ilev, self%state%lchnk) + wtp = carma_get_wght_pct(icol,ilev,self%state) end function wgtpct diff --git a/src/physics/carma/cam/carma_intr.F90 b/src/physics/carma/cam/carma_intr.F90 index 7315cbb70d..64c1d9942d 100644 --- a/src/physics/carma/cam/carma_intr.F90 +++ b/src/physics/carma/cam/carma_intr.F90 @@ -38,6 +38,8 @@ module carma_intr pbuf_get_index, pbuf_get_field, dtype_r8, pbuf_set_field use pio, only: var_desc_t + use wv_sat_methods, only: wv_sat_qsat_water + implicit none private @@ -89,6 +91,7 @@ module carma_intr public carma_get_wet_radius public carma_get_bin_rmass public carma_set_bin + public :: carma_get_wght_pct ! 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 @@ -202,9 +205,6 @@ module carma_intr ! elements. real (r8) :: carma_massscalefactor(NGROUP, NBIN) - ! sulfate weight percent -- updated in carma_timestep_tend - real(r8), allocatable, public, protected :: carma_sulfate_wght_pct(:,:,:) - contains @@ -807,14 +807,6 @@ subroutine carma_init(pbuf2d) call CARMAMODEL_InitializeModel(carma, lq_carma, pbuf2d, rc) if (rc < 0) call endrun('carma_init::CARMA_InitializeModel failed.') - ! allocate sulfate weight percent array - allocate(carma_sulfate_wght_pct(pcols,pver,begchunk:endchunk),stat=astat) - if( astat /= 0 ) then - write(iulog,*) 'carma_init: failed to allocate carma_sulfate_wght_pct, error = ',astat - call endrun - end if - carma_sulfate_wght_pct(:,:,:) = -huge(1._r8) - return end subroutine carma_init @@ -864,8 +856,6 @@ subroutine carma_final call CARMA_Destroy(carma, rc) if (rc < 0) call endrun('carma_final::CARMA_Destroy failed.') - deallocate(carma_sulfate_wght_pct) - return end subroutine carma_final @@ -1536,9 +1526,6 @@ subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rli ! Output diagnostic fields. call carma_output_diagnostics(state_loc, ptend, pbuf, cam_in, gpdiags, sbdiags, gsdiags, spdiags, bndiags) - ! save sulfate weight percent - carma_sulfate_wght_pct(:state%ncol,:,state%lchnk) = gsdiags(:state%ncol,:,I_GAS_H2SO4,GSDIAGS_WT) - end subroutine carma_timestep_tend !! Check the CARMA aerosol to make sure that for each aerosol the @@ -3837,6 +3824,8 @@ subroutine carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc) totvol(:ncol, :) = 0._r8 totmmr(:ncol, :) = 0._r8 + rhopdry(:ncol, :)= 0._r8 + rdry(:ncol, :) = 0._r8 do i = 1, nelems ielem = ielems(i) @@ -3852,12 +3841,16 @@ subroutine carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc) end do ! Add checks for totvol = 0 and nmr = 0 - rhopdry(:ncol, :) = totmmr(:ncol, :) / totvol(:ncol, :) + where(totvol(:ncol, :)>0._r8) + rhopdry(:ncol, :) = totmmr(:ncol, :) / totvol(:ncol, :) + end where 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) + where(nmr(:ncol, :)>0._r8) + rdry(:ncol, :) = ((3._r8 * totvol(:ncol, :) / nmr(:ncol, :)) / (4._r8 * PI)) ** (1._r8 / 3._r8) + end where return end subroutine carma_get_dry_radius @@ -4262,7 +4255,6 @@ end subroutine carma_get_total_mmr_cld !! @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 @@ -4297,13 +4289,19 @@ subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) if ((igroup < 1) .or. (igroup .gt. NGROUP)) then write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup rc = RC_ERROR - return + !return + end if + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR1: rc = ',rc) 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 + !return + end if + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR2: rc = ',rc) end if ! Get the constiuent index for water vapor (Q) @@ -4314,13 +4312,19 @@ subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) ! 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 + !if (rc < 0) return + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR3: rc = ',rc) + end if ! 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 + !if (rc < 0) return + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR4: rc = ',rc) + end if do icol = 1, ncol do iz = 1, pver @@ -4337,12 +4341,19 @@ subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) 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 - + if (rdry(icol, iz)>0._r8 .and. rhopdry(icol, iz)>0._r8) then + 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 + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR5: rc = ',rc) ! <====== + end if + else + rhopwet(icol,iz) = 0._r8 + rwet(icol, iz) = 0._r8 + end if else if (irhswell == I_PETTERS) then call carma_get_kappa(state, igroup, ibin, kappa, rc) if (rc < 0) return @@ -4351,12 +4362,19 @@ subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) 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 + !if (rc < 0) return + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR6: rc = ',rc) + end if 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 + !if (rc < 0) return + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR7: rc = ',rc) + end if + end if end do end do @@ -4365,6 +4383,10 @@ subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc) rwet(:ncol,:) = rwet(:ncol,:) / 1.e4_r8 rhopwet(:ncol,:) = rhopwet(:ncol,:) / 1.e3_r8 * 1.e6_r8 + if (rc/=RC_OK) then + call endrun('carma_get_wet_radius ERROR8: rc = ',rc) + end if + return end subroutine @@ -4437,4 +4459,40 @@ subroutine carma_get_bin_rmass(igroup, ibin, mass, rc) end subroutine carma_get_bin_rmass + function carma_get_wght_pct(icol,ilev,state) result(wtpct) + use sulfate_utils, only: wtpct_tabaz + + integer, intent(in) :: icol,ilev + type(physics_state), intent(in) :: state !! Physics state variables - before CARMA + + real(r8) :: wtpct + + integer :: rc !! return code + + real(r8) :: pvapl, es, qs, gc_cgs, rhoa + + rc = RC_OK + + ! Get relative humidity and vapor pressure + + call wv_sat_qsat_water(state%t(icol,ilev), state%pmid(icol,ilev), es, qs) + + ! seems like this is used in rwet calcuation to get pvapl water vapor pressure + + pvapl = es * 10._r8 / 1.e4_r8 + + ! need water vapor and gc_ptr is used above + + rhoa = (state%pmid(icol,ilev) * 10._r8) / (R_AIR * state%t(icol,ilev)) / 1.e3_r8 * 1.e6_r8 ! units ??? + + gc_cgs = state%q(icol,ilev,icnst4gas(carma%f_igash2o) ) * rhoa ! grams/cm3 ??? + + wtpct = wtpct_tabaz(carma, state%t(icol,ilev), gc_cgs, pvapl, rc) + + if (rc/=RC_OK) then + call endrun('carma_get_wght_pct: rc = ',rc) + end if + + end function carma_get_wght_pct + end module carma_intr