Merge pull request #266 from GEOS-ESM/feature/pcolarco/carbonChemLossUpd

Feature/pcolarco/carbon chem loss upd

CHANGELOG.md

@@ -16,6 +16,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Modified the file paths in carbon, sulfate, and nitrate ExtData.yaml files to used the revised version of the CEDS anthropogenic emissions. Note the previous version has an incorrect seasonal cycle.
 - Sulfate surface area density calculation in SU_Compute_Diags was incorrectly being passed the effective radius used for settling along with the sigma width of the number distribution. Properly it should be passed the number median radius, also present in the RC file. Added a hook to read that field from the RC file ("particle_radius_number"), store in SU grid comp, and pass to SU_Compute_Diags. This change is zero-diff to the SU internal state. It changes value of export SO4AREA.
 - Changed DMS concentration data holder from ExtData provided (SU_DMSO) to local copy (dmso_conc). This is relevant since if we run source tagged instances where we don't want DMS emissions we would zero out dmso_conc and that is what should be passed to DMSemission subroutine. This is zero diff except in that case.
+- It changes the formulation of the hydrophobic to hydrophilic conversion for carbon species,
+now defined by a time scale specified in the instance RC file. This is now specified by providing
+an e-folding time in days. This moves the time constant from outside the fortran to the run-time
+configurable RC file. This is not quite zero-diff with original code because of the precision of
+the specification, but testing shows nearly zero-diff result.
+
+- Also now present in the carbon instance RC files is a run-time configurable optional
+parameterized loss rate (e-folding time in days) per species and per mode. Default value for all is
+set to "-1" which means no use of this function.
 
 ### Fixed
 
@@ -32,7 +41,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Additional tuning parameters for the soil moisture and drylimit calculations for application specific tuning.
 - Required attributes for the 2D GOCART export fields in AERO_DP bundle have been set in subroutine append_to_bundle in Chem_AeroGeneric.F90. These export fields are imported by OBIO via Surface GC, and the missing of the attributes was causing the writing of surface import checkpoint to fail. The issue has been explained in detail on https://github.com/GEOS-ESM/GOCART/issues/258
-
 - Added export line to GOCART2G_GridCompMod to couple allow use of GOCART
   SU sulfate production tendency elsewhere in Chemistry, specifically for
   CARMA

ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90

@@ -55,6 +55,8 @@ module CA2G_GridCompMod
        real               :: fMonoterpenes = 0.0 ! Fraction of monoterpene emissions -> aerosol
        real               :: fIsoprene = 0.0 ! Franction of isoprene emissions -> aerosol
        real               :: fHydrophobic ! Initially hydrophobic portion
+       real               :: tConvPhobicToPhilic ! e-folding time [days] hydrophobic to hydrophilic
+       real               :: tChemLoss(2)        ! e-folding time [days] for parameterized chemistry loss
        logical            :: diurnal_bb   ! diurnal biomass burning
        real               :: eAircraftfuel       ! Aircraft emission factor: go from kg fuel to kg C
        real               :: aviation_layers(4)  ! heights of the LTO, CDS and CRS layers
@@ -163,6 +165,12 @@ subroutine SetServices (GC, RC)
 !   ----------------------------------------------
     call ESMF_ConfigGetAttribute (cfg, self%myDOW, label='my_day_of_week:', default=-1, __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%fhydrophobic, label='hydrophobic_fraction:', __RC__)
+    call ESMF_ConfigGetAttribute (cfg, self%tConvPhobicToPhilic, &
+                                  label='time_days_hydrophobic_to_hydrophilic:', default=2.5, __RC__)
+    call ESMF_ConfigFindLabel (cfg, 'time_days_chemical_destruction:', __RC__)
+    do i=1,size(self%tChemLoss)
+     call ESMF_ConfigGetAttribute (cfg, self%tChemLoss(i), default=-1., __RC__)
+    end do
     call ESMF_ConfigGetAttribute (cfg, self%ratPOM, label='pom_ca_ratio:', default=1.0, __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%fMonoterpenes, label='monoterpenes_emission_fraction:', default=0.0, __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%fIsoprene, label='isoprene_emission_fraction:', default=0.0, __RC__)
@@ -998,9 +1006,20 @@ subroutine Run2 (GC, import, export, clock, RC)
     end if
 
 !   Ad Hoc transfer of hydrophobic to hydrophilic aerosols
-!   Following Chin's parameterization, the rate constant is
-!   k = 4.63e-6 s-1 (.4 day-1; e-folding time = 2.5 days)
-    call phobicTophilic (intPtr_phobic, intPtr_philic, HYPHIL, self%km, self%cdt, MAPL_GRAV, delp, __RC__)
+!   Rate controlled in RC file; tConvPhobicToPhilic < 0 means no transfer
+    call phobicToPhilic (intPtr_phobic, intPtr_philic, HYPHIL, &
+                         self%tConvPhobicToPhilic, self%km, self%cdt, MAPL_GRAV, delp, __RC__)
+
+!   Ad Hoc chemical destruction of carbon
+!   This applies a simple exponential decay to both hydrophobic and
+!   hydrophilic modes with the time constant tChemLoss (e-folding
+!   time in days)
+    do n = 1, self%nbins
+       call MAPL_VarSpecGet(InternalSpec(n), SHORT_NAME=short_name, __RC__)
+       call MAPL_GetPointer(internal, NAME=short_name, ptr=int_ptr, __RC__)
+       call carbonChemLoss (self%km, self%klid, n, self%cdt, MAPL_GRAV, delp, &
+                            self%tChemLoss(n), int_ptr, CH, __RC__)
+    end do
 
 !   CA Settling
 !   -----------

ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_StateSpecs.rc

@@ -89,6 +89,7 @@ category: EXPORT
  *DP          | kg m-2 s-1   | xy    | N     | nbins                          | * Aerosol Dry Deposition (Bin %d)
  *WT          | kg m-2 s-1   | xy    | N     | nbins                          | * Aerosol Wet Deposition (Bin %d)
  *SV          | kg m-2 s-1   | xy    | N     | nbins                          | * Aerosol Convective Scavenging (Bin %d)
+ *CH          | kg m-2 s-1   | xy    | N     | nbins                          | * Aerosol Parameterized Chemical Loss (Bin %d)
  *EMAN        | kg m-2 s-1   | xy    | N     |                                | * Aerosol Anthropogenic Emissions
  *EMBB        | kg m-2 s-1   | xy    | N     |                                | * Aerosol Biomass Burning Emissions
  *EMBF        | kg m-2 s-1   | xy    | N     |                                | * Aerosol Biofuel Emissions

ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_instance_CA.bc.rc

@@ -14,6 +14,12 @@ aviation_vertical_layers: 0.0 100.0 9.0e3 10.0e3
 # Initially hydrophobic portion
 hydrophobic_fraction: 0.8
 
+# Rate of conversion of hydrophobic to hydrophilic [days]
+time_days_hydrophobic_to_hydrophilic: 2.5
+
+# Rate of chemical destruction of carbon species [days]
+time_days_chemical_destruction: -1.  -1.
+
 # Scavenging efficiency per bin [km-1] (NOT USED UNLESS RAS IS CALLED)
 fscav: 0.0  0.4 
 

ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_instance_CA.br.rc

@@ -20,6 +20,12 @@ particle_radius_microns: 0.35 0.35
 # Initially hydrophobic portion
 hydrophobic_fraction: 0.5
 
+# Rate of conversion of hydrophobic to hydrophilic [days]
+time_days_hydrophobic_to_hydrophilic: 2.5
+
+# Rate of chemical destruction of carbon species [days]
+time_days_chemical_destruction: -1.   -1.
+
 # Scavenging efficiency per bin [km-1] (NOT USED UNLESS RAS IS CALLED)
 fscav: 0.0  0.4 
 

ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_instance_CA.oc.rc

@@ -24,6 +24,12 @@ particle_radius_microns: 0.35 0.35
 # Initially hydrophobic portion
 hydrophobic_fraction: 0.5
 
+# Rate of conversion of hydrophobic to hydrophilic [days]
+time_days_hydrophobic_to_hydrophilic: 2.5
+
+# Rate of chemical destruction of carbon species [days]
+time_days_chemical_destruction: -1.  -1.
+
 # Scavenging efficiency per bin [km-1] (NOT USED UNLESS RAS IS CALLED)
 fscav: 0.0  0.4 
 

Process_Library/GOCART2G_Process.F90

@@ -41,7 +41,8 @@ module  GOCART2G_Process
    public wetRadius
    public hoppelCorrection
    public CAEmission
-   public phobicTophilic
+   public phobicToPhilic
+   public carbonChemLoss
    public NIheterogenousChem
    public SulfateDistributeEmissions
    public DMSemission
@@ -4682,16 +4683,17 @@ end subroutine distribute_aviation_emissions
 
 !BOP
 !
-! !IROUTINE: phobicTophilic
+! !IROUTINE: phobicToPhilic
 !
 ! !INTERFACE:
    subroutine phobicTophilic (aerosol_phobic, aerosol_philic, aerosol_toHydrophilic, &
-                              km, cdt, grav, delp, rc)
+                              tConvPhobicToPhilic, km, cdt, grav, delp, rc)
 
 ! !USES:
    implicit NONE
 
 ! !INPUT PARAMETERS:
+   real, intent(in)      :: tConvPhobicToPhilic  ! e-folding time in days to transfer
    integer, intent(in)   :: km   ! total model level
    real, intent(in)      :: cdt  ! chemistry model time-step [sec]
    real, intent(in)      :: grav ! [m/sec^2]
@@ -4706,18 +4708,25 @@ subroutine phobicTophilic (aerosol_phobic, aerosol_philic, aerosol_toHydrophilic
 
 ! !Local Variables
    integer :: i, j, k
-   real :: qUpdate, delq
+   real :: qUpdate, delq, ts
 
 !EOP
 !------------------------------------------------------------------------------------
 !  Begin...
 
    if(associated(aerosol_toHydrophilic)) aerosol_toHydrophilic = 0.0
 
+!  tConvPhobicToPhilic is the e-folding time (in days) of the conversion
+!  If < 0 no conversion is desired; exit the subroutine
+   if(tConvPhobicToPhilic < 0) then
+    __RETURN__(__SUCCESS__)
+   endif
+   ts = tConvPhobicToPhilic*86400.
+
    do k = 1, km
     do j = 1, ubound(delp, 2)
      do i = 1, ubound(delp, 1)
-      qUpdate = aerosol_phobic(i,j,k)*exp(-4.63e-6*cdt)
+      qUpdate = aerosol_phobic(i,j,k)*exp(-cdt/ts)
       qUpdate = max(qUpdate,1.e-32)
       delq = max(0.,aerosol_phobic(i,j,k)-qUpdate)
       aerosol_phobic(i,j,k) = qUpdate
@@ -4732,6 +4741,70 @@ subroutine phobicTophilic (aerosol_phobic, aerosol_philic, aerosol_toHydrophilic
    __RETURN__(__SUCCESS__)
   end subroutine phobicTophilic
 
+!============================================================================
+!BOP
+!
+! !IROUTINE: carbonChemLoss
+!
+! !INTERFACE:
+   subroutine carbonChemLoss (km, klid, n, cdt, grav, delp, &
+                              tChemLoss, int_qa, fluxout, rc)
+
+! !USES:
+   implicit NONE
+
+! !INPUT PARAMETERS:
+   real, intent(in)       :: tChemLoss  ! e-folding loss time [days]
+   integer, intent(in)    :: km         ! total model levels
+   integer, intent(in)    :: klid       ! index for pressure lid
+   integer, intent(in)    :: n          ! bin index number
+   real, intent(in)       :: cdt        ! time step [s]
+   real, intent(in)       :: grav       ! acceleration of gravity [m/sec^2]
+   real, dimension(:,:,:), intent(inout)          :: int_qa  ! aerosol [kg/kg]
+   real, pointer, dimension(:,:,:), intent(in)    :: delp    ! pressure level thickness [Pa]
+
+! !OUTPUT PARAMETERS:
+
+   real, pointer, dimension(:,:,:), intent(inout) :: fluxout ! Mass lost by chemistry [kg/m^2/s]
+   integer, optional, intent(out)                 :: rc      ! Error return code:
+                                                             !  0 - all is well
+                                                             !  1 -
+
+! !Local Variables
+   integer :: i, j, k
+   real :: qUpdate, delq, ts
+
+!EOP
+!------------------------------------------------------------------------------------
+!  Begin...
+
+   if(associated(fluxout)) fluxout(:,:,n) = 0.0
+
+!  tChemLoss is the e-folding time (in days) of parameterized chemistry loss
+!  If < 0 no loss is desired; exit the subroutine
+   if(tChemLoss < 0) then
+    __RETURN__(__SUCCESS__)
+   endif
+   ts = tChemLoss*86400.
+
+   do k = klid, km
+    do j = 1, ubound(delp, 2)
+     do i = 1, ubound(delp, 1)
+      qUpdate = int_qa(i,j,k)*exp(-cdt/ts)
+      qUpdate = max(qUpdate,1.e-32)
+      delq = max(0.,int_qa(i,j,k)-qUpdate)
+      int_qa(i,j,k) = qUpdate
+      if(associated(fluxout)) &
+       fluxout(i,j,n) = fluxout(i,j,n) &
+        + delq*delp(i,j,k)/grav/cdt
+
+     end do
+    end do
+   end do
+
+   __RETURN__(__SUCCESS__)
+  end subroutine carbonChemLoss
+
 
 !============================================================================
 !BOP