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
fvitt · Nov 14, 2023 · 44ed90f · 44ed90f
1 parent 326a6c3
commit 44ed90f
Show file tree

Hide file tree

Showing 3 changed files with 105 additions and 41 deletions.
diff --git 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
@@ -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
@@ -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