Merge pull request #249 from GEOS-ESM/develop to main

Updates to CARMA, and QFED paths

CARMAchem_GridComp/CARMA/source/base/carma_constants_mod.F90

@@ -47,9 +47,15 @@ module carma_constants_mod
 
 	!! Define molecular weight of sulphuric acid [ g / mole ]
 	real(kind=f), parameter :: WTMOL_H2SO4 = 98.078479_f
+
+        !! Define molecular weight of sulfur dioxide [ g / mole ]
+	real(kind=f), parameter :: WTMOL_SO2 = 64.066_f
 
-  !! Define reference pressure, e.g. for potential temp calcs [ dyne / cm^2 ]
-  real(kind=f), parameter :: PREF = 1000.e+3_f
+        !! Define molecular weight of nitric acid [ g / mole ]
+	real(kind=f), parameter :: WTMOL_HNO3 = 62.996_f
+
+        !! Define reference pressure, e.g. for potential temp calcs [ dyne / cm^2 ]
+        real(kind=f), parameter :: PREF = 1000.e+3_f
 
 	!! Define conversion factor for mb to cgs [ dyne / cm^2 ] units
 	real(kind=f), parameter :: RMB2CGS = 1000.e+0_f

CARMAchem_GridComp/CARMA/source/base/carma_enums_mod.F90

@@ -65,6 +65,7 @@ module carma_enums_mod
   integer, public, parameter :: I_FITZGERALD   = 1   !! Fitzgerald
   integer, public, parameter :: I_GERBER       = 2   !! Gerber
   integer, public, parameter :: I_WTPCT_H2SO4  = 3   !! The weight percent method for sulfate aerosol
+  integer, public, parameter :: I_WTPCT_STS    = 4   !! The weight percent method for sts
 
   !  Define vallues of flag used for particle swelling composition (Fiztgerald)
   integer, public, parameter :: I_SWF_NH42SO4   = 1   !! (NH4)2SO4
@@ -88,6 +89,7 @@ module carma_enums_mod
   integer, public, parameter :: I_VAPRTN_H2O_MURPHY2005    = 2   !! H2O, Murphy & Koop [2005]
   integer, public, parameter :: I_VAPRTN_H2O_GOFF1946      = 3   !! H2O, Goff & Gratch [1946], used in CAM
   integer, public, parameter :: I_VAPRTN_H2SO4_AYERS1980   = 4   !! H2SO4, Ayers [1980] & Kumala [1990]
+  integer, public, parameter :: I_VAPRTN_NULL              = 5   !! For non-condensing gases
 
   ! Routines to calculate fall velocities
   integer, public, parameter :: I_FALLRTN_STD              = 1   !! Standard CARMA 2.3 routine (spherical only)
@@ -103,6 +105,7 @@ module carma_enums_mod
   integer, public, parameter :: I_GCOMP_H2O             = 1   !! Water Vapor
   integer, public, parameter :: I_GCOMP_H2SO4           = 2   !! Sulphuric Acid
   integer, public, parameter :: I_GCOMP_SO2             = 3   !! Sulfer Dioxide
+  integer, public, parameter :: I_GCOMP_HNO3            = 4   !! Nitric Acid
 
   ! How is the CARMA group represented in the parent model
   integer, public, parameter :: I_CNSTTYPE_PROGNOSTIC   = 1   !! Prognostic, advected constituent for each bin

CARMAchem_GridComp/CARMA/source/base/carma_globaer.h

@@ -180,6 +180,7 @@
 #define dt_threshold  carma%f_dt_threshold
 #define igash2o       carma%f_igash2o
 #define igash2so4     carma%f_igash2so4
+#define igashno3      carma%f_igashno3
 #define igasso2       carma%f_igasso2
 #define tstick        carma%f_tstick
 #define gsticki       carma%f_gsticki

CARMAchem_GridComp/CARMA/source/base/carma_mod.F90

@@ -850,6 +850,7 @@ subroutine CARMA_InitializeGrowth(carma, rc)
     carma%f_igash2o   = 0
     carma%f_igash2so4 = 0
     carma%f_igasso2   = 0
+    carma%f_igashno3  = 0
 
     do igas = 1, carma%f_NGAS
       if (carma%f_gas(igas)%f_icomposition == I_GCOMP_H2O) then
@@ -858,6 +859,8 @@ subroutine CARMA_InitializeGrowth(carma, rc)
         carma%f_igash2so4 = igas
       else if (carma%f_gas(igas)%f_icomposition == I_GCOMP_SO2) then
         carma%f_igasso2 = igas
+      else if (carma%f_gas(igas)%f_icomposition == I_GCOMP_HNO3) then
+        carma%f_igashno3 = igas
       end if
     end do
 

CARMAchem_GridComp/CARMA/source/base/carma_types_mod.F90

@@ -401,6 +401,7 @@ module carma_types_mod
     integer                                       :: f_igash2o
     integer                                       :: f_igash2so4
     integer                                       :: f_igasso2
+    integer                                       :: f_igashno3
     integer                                       :: f_maxsubsteps 
     integer                                       :: f_minsubsteps 
     integer                                       :: f_maxretries 
@@ -817,4 +818,4 @@ module carma_types_mod
     real(kind=f), allocatable, dimension(:,:,:)  :: f_dtpart     ! (NZ,NBIN,NGROUP)
     real(kind=f)                                 :: f_phprod
    end type carmastate_type   
-end module
+end module

CARMAchem_GridComp/CARMA/source/base/nsubsteps.F90

@@ -94,15 +94,18 @@ subroutine nsubsteps(carma, cstate, iz, dtime_save, ntsubsteps, rc)
           elseif( itype(iepart) .eq. I_VOLATILE ) then
 
             do ibin = NBIN-1, 1, -1
-!              if( pc(iz,ibin,iepart) / xmet(iz) / ymet(iz) / zmet(iz) .gt. conmax * pconmax(iz,ig) )then
-!write(*,'(3(i3,2x),6(e15.7,2x))') iz, ibin, iepart, pc(iz,ibin,iepart), conmax,  pconmax(iz,ig), xmet(iz), ymet(iz), zmet(iz)
-!              if( pc(iz,ibin,iepart) .gt. conmax * pconmax(iz,ig) * xmet(iz) * ymet(iz) * zmet(iz) )then
-! Check the presence of NaN or Inf
-  if(ieee_is_nan(pc(iz,ibin,iepart))) print *, '     pc isnan', iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
-  if(ieee_is_nan(pconmax(iz,ig)))    print *,  'pconmax isnan', iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
-  if(.not. ieee_is_finite(pc(iz,ibin,iepart))) print *, '     pc isinf', iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
-  if(.not. ieee_is_finite(pconmax(iz,ig)))     print *, 'pconmax isinf', iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
-              if( pc(iz,ibin,iepart)/pconmax(iz,ig) .gt. conmax * xmet(iz) * ymet(iz) * zmet(iz) )then
+              ! Check the presence of NaN or Inf
+              ! PAC: Why is this done with "print" instead of "write" like
+              !  the rest of the error messages?
+              if(ieee_is_nan(pc(iz,ibin,iepart))) print *, 'pc isnan', iz, &
+                ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
+              if(ieee_is_nan(pconmax(iz,ig)))    print *,  'pconmax isnan', iz, &
+                ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
+              if(.not. ieee_is_finite(pc(iz,ibin,iepart))) print *, 'pc isinf', &
+                iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
+              if(.not. ieee_is_finite(pconmax(iz,ig))) print *, 'pconmax isinf', &
+                iz, ibin, iepart, pc(iz,ibin,iepart), pconmax(iz,ig)
+              if( pc(iz,ibin,iepart)/pconmax(iz,ig) .gt. conmax * xmet(iz) * ymet(iz) * zmet(iz))then
                 ibin_small(ig) = ibin
               endif
             enddo

CARMAchem_GridComp/CARMA/source/base/rhopart.F90

@@ -45,6 +45,9 @@ subroutine rhopart(carma, cstate, rc)
   real(kind=f)    :: mcore(NBIN)
   real(kind=f)    :: r_ratio
   real(kind=f)    :: h2o_mass
+  real(kind=f)    :: h2o_vmr
+  real(kind=f)    :: hno3_vmr
+  real(kind=f)    :: h2so4m
 
   1 format(/,'rhopart::WARNING - core mass > total mass, truncating : iz=',i4,',igroup=',&
               i4,',ibin=',i4,',total mass=',e10.3,',core mass=',e10.3,',using rhop=',f9.4)
@@ -118,16 +121,25 @@ subroutine rhopart(carma, cstate, rc)
       ! and dry radius are the same.
 
       ! Determine the weight percent of sulfate, and store it for later use.
-      if (irhswell(igroup) == I_WTPCT_H2SO4) then
+      if (irhswell(igroup) == I_WTPCT_H2SO4 .OR. &
+          (irhswell(igroup) == I_WTPCT_STS .AND. t(iz) > 200._f)) then
         h2o_mass     = gc(iz, igash2o) / (xmet(iz) * ymet(iz) * zmet(iz))
+      else if (irhswell(igroup) == I_WTPCT_STS) then
+        ! h2o_vmr and hno3_vmr are taken from gas concentrations while h2so4m
+        !   needs to be the mass of h2so4 in particles and in gas phase
+        h2o_vmr     = gc(iz, igash2o) / gwtmol(igash2o) * WTMOL_AIR
+        hno3_vmr    = gc(iz, igashno3) / gwtmol(igashno3) * WTMOL_AIR 
+        h2so4m       = sum(rmass(:, igroup) * pc(iz, :, iepart)) + &
+                         gc(iz, igash2so4) / (xmet(iz) * ymet(iz) * zmet(iz))
       end if
 
       ! Loop over particle size bins.
       do ibin = 1,NBIN
 
         ! If humidity affects the particle, then determine the equilbirium
         ! radius and density based upon the relative humidity.
-        if (irhswell(igroup) == I_WTPCT_H2SO4) then
+        if (irhswell(igroup) == I_WTPCT_H2SO4 .OR. &
+            (irhswell(igroup) == I_WTPCT_STS .AND. t(iz) > 200._f)) then
 
           ! rlow
           call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &
@@ -147,6 +159,31 @@ subroutine rhopart(carma, cstate, rc)
             h2o_vp=pvapl(iz, igash2o), temp=t(iz))
           if (rc < 0) return
 
+        ! If STS parameterization is selected, use hno3, h2o, and h2so4 mass
+        ! to determine equilibrium radius and density
+        else if (irhswell(igroup) == I_WTPCT_STS) then
+
+          ! rlow
+          call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &
+            rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_vmr=h2o_vmr, &
+            h2o_vp=pvapl(iz, igash2o), temp=t(iz), press=cstate%f_p(iz), h2so4m = h2so4m, &
+            hno3_vmr = hno3_vmr)
+          if (rc < 0) return
+
+          ! rup
+          call getwetr(carma, igroup, relhum(iz), rup(ibin,igroup), rup_wet(iz,ibin,igroup), &
+            rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_vmr=h2o_vmr, &
+            h2o_vp=pvapl(iz, igash2o), temp=t(iz), press=cstate%f_p(iz), h2so4m = h2so4m, &
+            hno3_vmr = hno3_vmr)
+          if (rc < 0) return
+
+          ! r
+          call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
+            rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_vmr=h2o_vmr, &
+            h2o_vp=pvapl(iz, igash2o), temp=t(iz), press=cstate%f_p(iz), h2so4m = h2so4m, &
+            hno3_vmr = hno3_vmr)
+          if (rc < 0) return
+
         else
           ! rlow
           call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &

CARMAchem_GridComp/CARMA/source/base/setupgkern.F90

@@ -124,6 +124,15 @@ subroutine setupgkern(carma, cstate, rc)
         ! for water vapor.
         akelvini(k, igas) = akelvini(k, igash2o)
       end do   
+    else if (igas .eq. igashno3) then
+      ! Not growing anything with HNO3
+      ! Just grabbing akelvin for h2so4(liq)/water vapor(ice)
+      do k = 1, NZ  
+        surf_tens = sulfate_surf_tens(carma, wtpct(k), t(k), rc)
+        rho_H2SO4 = sulfate_density(carma, wtpct(k), t(k), rc)
+        akelvin(k, igas) = 2._f * gwtmol(igas) * surf_tens / (t(k) * rho_H2SO4 * RGAS)
+        akelvini(k, igas) = akelvini(k, igash2o)
+      end do
     else 
 
       ! Condensing gas is not yet configured.

CARMAchem_GridComp/CARMA/source/base/setupgrow.F90

@@ -99,6 +99,19 @@ subroutine setupgrow(carma, cstate, rc)
       rlhe(k,igash2so4) = rlhe(k, igash2o)
       rlhm(k,igash2so4) = rlhe(k, igash2o)
     end if
+
+    ! Properties for HNO3
+    ! PAC: Making these the same as H2SO4 for now.
+    if (igashno3 /= 0) then
+      ! Diffusivity
+      rhoa_cgs = rhoa(k) / (xmet(k) * ymet(k) * zmet(k))
+      aden     = rhoa_cgs * AVG / WTMOL_AIR
+      diffus(k,igashno3) = 1.76575e+17_f * sqrt(t(k)) / aden
+
+      ! HACK: make H2SO4 latent heats same as water
+      rlhe(k,igashno3) = rlhe(k, igash2o)
+      rlhm(k,igashno3) = rlhe(k, igash2o)
+    end if
 
   enddo
 

CARMAchem_GridComp/CARMA/source/base/setupvf.F90

@@ -29,8 +29,10 @@ subroutine setupvf(carma, cstate, rc)
 
   ! Local declarations  
   integer                 :: igroup, i, j, k, k1, k2, ibin, iz, nzm1
+  integer                 :: iepart
   real(kind=f)            :: r_tmp(NBIN,NGROUP)
   real(kind=f)            :: h2o_mass
+  real(kind=f)            :: h2o_vmr, hno3_vmr, h2so4m
 
   !  Define formats
   2 format(/,'Fall velocities and Reynolds number (prior to interpolation)')
@@ -42,19 +44,31 @@ subroutine setupvf(carma, cstate, rc)
   do igroup = 1, NGROUP
 
     ! Special handling for vf_const < 0 (if, then abs(vf_const) = dry particle radius [cm])
+    ! PAC: do I need to add logic for STS here?
     if( ifall == 0 .and. vf_const < 0._f) then
       r_tmp = abs(vf_const)
       do iz = 1, NZ
        do ibin = 1, NBIN
-        if (irhswell(igroup) == I_WTPCT_H2SO4) then
+         if (irhswell(igroup) == I_WTPCT_H2SO4 .OR. &
+             (irhswell(igroup) == I_WTPCT_STS .AND. t(iz) > 200.)) then
           h2o_mass = gc(iz, igash2o) / (xmet(iz) * ymet(iz) * zmet(iz))
-          call getwetr(carma, igroup, relhum(iz), r_tmp(ibin,igroup), r_wet(iz,ibin,igroup), &
+          call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
                        rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
                        h2o_vp=pvapl(iz, igash2o), temp=t(iz))
-        else
-          call getwetr(carma, igroup, relhum(iz), r_tmp(ibin,igroup), r_wet(iz,ibin,igroup), &
+         else if (irhswell(igroup) == I_WTPCT_STS) then
+          iepart = ienconc(igroup)     ! element of particle number concentration 
+          h2o_vmr     = gc(iz, igash2o) / gwtmol(igash2o) * WTMOL_AIR
+          hno3_vmr    = gc(iz, igashno3) / gwtmol(igashno3) * WTMOL_AIR 
+          h2so4m       = sum(rmass(:, igroup) * pc(iz, :, iepart)) + &
+                             gc(iz, igash2so4) / (xmet(iz) * ymet(iz) * zmet(iz))
+          call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
+            rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_vmr=h2o_vmr, &
+            h2o_vp=pvapl(iz, igash2o), temp=t(iz), press=cstate%f_p(iz), h2so4m = h2so4m, &
+            hno3_vmr = hno3_vmr)
+         else
+          call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
                        rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc)
-        endif
+         endif
         if (rc < 0) return
        enddo
       enddo
@@ -90,11 +104,22 @@ subroutine setupvf(carma, cstate, rc)
       ! Special handling for vf_const < 0 (if, then abs(vf_const) = dry particle radius [cm])
       do iz = 1, NZ
        do ibin = 1, NBIN
-         if (irhswell(igroup) == I_WTPCT_H2SO4) then
+         if (irhswell(igroup) == I_WTPCT_H2SO4 .OR. &
+             (irhswell(igroup) == I_WTPCT_STS .AND. t(iz) > 200.)) then
           h2o_mass = gc(iz, igash2o) / (xmet(iz) * ymet(iz) * zmet(iz))
           call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
                        rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
                        h2o_vp=pvapl(iz, igash2o), temp=t(iz))
+         else if (irhswell(igroup) == I_WTPCT_STS) then
+          iepart = ienconc(igroup)     ! element of particle number concentration 
+          h2o_vmr     = gc(iz, igash2o) / gwtmol(igash2o) * WTMOL_AIR
+          hno3_vmr    = gc(iz, igashno3) / gwtmol(igashno3) * WTMOL_AIR 
+          h2so4m       = sum(rmass(:, igroup) * pc(iz, :, iepart)) + &
+                             gc(iz, igash2so4) / (xmet(iz) * ymet(iz) * zmet(iz))
+          call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
+            rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_vmr=h2o_vmr, &
+            h2o_vp=pvapl(iz, igash2o), temp=t(iz), press=cstate%f_p(iz), h2so4m = h2so4m, &
+            hno3_vmr = hno3_vmr)
          else
           call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
                        rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc)