prc: updates to include separate explosive and degassing volcanic emi…

…ssions for SO2

CHANGELOG.md

@@ -13,6 +13,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+## [Unreleased - [email protected]] - 2024-01-08
+
+### Changed
+
+- Changed SU2G_instance_SU.rc to now have separate filename inputs for explosive and degassing volcanoes
+- Moved present volcanic emission inventories to one or the other line for these new entries; set other
+  line /dev/null; this is stop gap until next time we update volcanic emission inventories, at which 
+  point will provide (for AMIP and AMIP.20C) separate explosive and degassing emissions
+- Made accommodating changes for above in SU2G_GridCompMod.F90 and in the Process Library
+- Verified zero diff in current configuration (this is true of tracers and restarts, but not diagnostics:
+  until an actual split is made in the input emissions then the volcanic emissions are being assigned to
+  one or the other emission diagnostics (explosive or degassing).
+
 ## [v2.2.1] - 2023-05-30
 
 ### Fixed

ESMF/GOCART2G_GridComp/SU2G_GridComp/AMIP.20C/SU2G_instance_SU.rc

@@ -8,7 +8,8 @@ aerosol_monochromatic_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opti
 nbins: 4
 
 # Volcanic pointwise sources
-volcano_srcfilen: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_Carns.%y4%m2%d2.rc
+volcano_srcfilen_explosive: /dev/null
+volcano_srcfilen_degassing: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_Carns.%y4%m2%d2.rc
 
 # Heights [m] of LTO, CDS and CRS aviation emissions layers
 aviation_vertical_layers: 0.0 100.0 9.0e3 10.0e3

ESMF/GOCART2G_GridComp/SU2G_GridComp/AMIP/SU2G_instance_SU.rc

@@ -8,7 +8,8 @@ aerosol_monochromatic_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opti
 nbins: 4
 
 # Volcanic pointwise sources
-volcano_srcfilen: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_Carns.%y4%m2%d2.rc
+volcano_srcfilen_explosive: /dev/null
+volcano_srcfilen_degassing: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_Carns.%y4%m2%d2.rc
 
 # Heights [m] of LTO, CDS and CRS aviation emissions layers
 aviation_vertical_layers: 0.0 100.0 9.0e3 10.0e3

ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90

@@ -54,6 +54,9 @@ module SU2G_GridCompMod
 !===========================================================================
 !  !Sulfer state
    type :: ThreadWorkspace
+      integer :: nymd_last = -1 ! Previous nymd. Updated daily
+
+!     Degassing volcanoes
       integer :: nVolc = 0
       real, allocatable, dimension(:)  :: vLat, &
                                           vLon, &
@@ -62,7 +65,17 @@ module SU2G_GridCompMod
                                           vCloud
       integer, allocatable, dimension(:) :: vStart, &
                                             vEnd
-      integer :: nymd_last = -1 ! Previous nymd. Updated daily
+!     Explosive volcanoes
+      integer :: nVolcE = 0
+      real, allocatable, dimension(:)  :: vLatE, &
+                                          vLonE, &
+                                          vSO2E, &
+                                          vElevE, &
+                                          vCloudE
+      integer, allocatable, dimension(:) :: vStartE, &
+                                            vEndE
+
+!     Other point emissions of Sulfate (SO4)
       integer                         :: nPts = -1
       integer, allocatable, dimension(:)  :: pstart, pend
       real, allocatable, dimension(:)     :: pLat, &
@@ -86,8 +99,10 @@ module SU2G_GridCompMod
       !real, pointer :: h2o2_init(:,:,:)
 
 !     Special handling for volcanic emissions
-      character(len=255) :: volcano_srcfilen
-!     !Workspae for point emissions
+      character(len=255) :: volcano_srcfilen_degassing
+      character(len=255) :: volcano_srcfilen_explosive
+
+!     Workspace for point emissions
       logical                :: doing_point_emissions = .false.
       character(len=255)     :: point_emissions_srcfilen   ! filename for pointwise emissions
       type(ThreadWorkspace), allocatable :: workspaces(:)
@@ -174,7 +189,8 @@ subroutine SetServices ( GC, RC )
     allocate(self%sigma(self%nbins), __STAT__)
 
 !   process SU-specific items
-    call ESMF_ConfigGetAttribute(cfg, self%volcano_srcfilen, label='volcano_srcfilen:', __RC__)
+    call ESMF_ConfigGetAttribute(cfg, self%volcano_srcfilen_degassing, label='volcano_srcfilen_degassing:', __RC__)
+    call ESMF_ConfigGetAttribute(cfg, self%volcano_srcfilen_explosive, label='volcano_srcfilen_explosive:', __RC__)
     call ESMF_ConfigGetAttribute(cfg, self%eAircraftFuel, label='aircraft_fuel_emission_factor:', __RC__)
     call ESMF_ConfigGetAttribute(cfg, self%fSO4anth, label='so4_anthropogenic_fraction:', __RC__)
     call ESMF_ConfigGetAttribute(cfg, self%sigma, label='sigma:', __RC__)
@@ -718,7 +734,7 @@ subroutine Run1 (GC, import, export, clock, RC)
     real, dimension(:,:), allocatable :: so2biomass_src, so2biomass_src_, so2anthro_l1_src, &
                                          so2anthro_l2_src, so2ship_src, so4ship_src, dmso_conc, &
                                          aviation_lto_src, aviation_cds_src, aviation_crs_src
-    integer, dimension(:), allocatable  :: iPointVolc, jPointVolc, iPoint, jPoint
+    integer, dimension(:), allocatable  :: iPoint, jPoint
     real, dimension(:,:,:), allocatable :: emissions_point
     character (len=ESMF_MAXSTR)  :: fname ! file name for point source emissions
     logical :: fileExists
@@ -812,34 +828,48 @@ subroutine Run1 (GC, import, export, clock, RC)
     where(1.01*aviation_cds_src > MAPL_UNDEF ) aviation_cds_src = 0.
     where(1.01*aviation_crs_src > MAPL_UNDEF ) aviation_crs_src = 0.
 
+!   Start with a clean emission diagnostic
+    if(associated(SUEM)) SUEM = 0.0    
+
+
 !   Update emissions/production if necessary (daily)
 !   -----------------------------------------------
     thread = MAPL_get_current_thread()
     workspace => self%workspaces(thread)
 
+!   Update Volcanic SO2 Emissions Daily
     if(workspace%nymd_last /= nymd) then
        workspace%nymd_last = nymd
 
-!      Get pointwise SO2 and altitude of volcanoes from a daily file data base
-       if(index(self%volcano_srcfilen,'volcanic_') /= 0) then
-          call StrTemplate(fname, self%volcano_srcfilen, xid='unknown', &
+!      DEGASSING: Get pointwise SO2 and altitude of volcanoes from a daily file data base
+       workspace%nVolc = 0  ! case of /dev/null (no volcanoes) or ill-formed filename
+       if(index(self%volcano_srcfilen_degassing,'volcanic_') /= 0) then
+          call StrTemplate(fname, self%volcano_srcfilen_degassing, xid='unknown', &
                             nymd=nymd, nhms=120000 )
           call ReadPointEmissions (nymd, fname, workspace%nVolc, workspace%vLat, workspace%vLon, &
                                    workspace%vElev, workspace%vCloud, workspace%vSO2, workspace%vStart, &
                                    workspace%vEnd, label='volcano', __RC__)
           workspace%vSO2 = workspace%vSO2 * fMassSO2 / fMassSulfur
-!         Special possible case
-          if(self%volcano_srcfilen(1:9) == '/dev/null') workspace%nVolc = 0
+       end if
+
+!      EXPLOSIVE: Get pointwise SO2 and altitude of volcanoes from a daily file data base
+       workspace%nVolcE = 0  ! case of /dev/null (no volcanoes) or ill-formed filename
+       if(index(self%volcano_srcfilen_explosive,'volcanic_') /= 0) then
+          call StrTemplate(fname, self%volcano_srcfilen_explosive, xid='unknown', &
+                            nymd=nymd, nhms=120000 )
+          call ReadPointEmissions (nymd, fname, workspace%nVolcE, workspace%vLatE, workspace%vLonE, &
+                                   workspace%vElevE, workspace%vCloudE, workspace%vSO2E, workspace%vStartE, &
+                                   workspace%vEndE, label='volcano', __RC__)
+          workspace%vSO2 = workspace%vSO2 * fMassSO2 / fMassSulfur
        end if
     end if
 
-!   Apply volcanic emissions
-!   ------------------------
+!   DEGASSING: Apply volcanic emissions
+!   -----------------------------------
     if (workspace%nVolc > 0) then
-       if (associated(SO2EMVE)) SO2EMVE=0.0
        if (associated(SO2EMVN)) SO2EMVN=0.0
-       allocate(iPointVolc(workspace%nVolc), jPointVolc(workspace%nVolc),  __STAT__)
-       call MAPL_GetHorzIJIndex(workspace%nVolc, iPointVolc, jPointVolc, &
+       allocate(iPoint(workspace%nVolc), jPoint(workspace%nVolc),  __STAT__)
+       call MAPL_GetHorzIJIndex(workspace%nVolc, iPoint, jPoint, &
                                 grid = grid,               &
                                 lon  = workspace%vLon/real(MAPL_RADIANS_TO_DEGREES), &
                                 lat  = workspace%vLat/real(MAPL_RADIANS_TO_DEGREES), &
@@ -850,8 +880,30 @@ subroutine Run1 (GC, import, export, clock, RC)
            end if
 
        call SUvolcanicEmissions (workspace%nVolc, workspace%vStart, workspace%vEnd, workspace%vSO2, workspace%vElev, &
-                                 workspace%vCloud, iPointVolc, jPointVolc, nhms, SO2EMVN, SO2EMVE, SO2, nSO2, SUEM, &
+                                 workspace%vCloud, iPoint, jPoint, nhms, SO2EMVN, SO2, nSO2, SUEM, &
                                  self%km, self%cdt, MAPL_GRAV, zle, delp, area, workspace%vLat, workspace%vLon, __RC__)
+       deallocate(iPoint, jPoint, __STAT__)
+    end if
+
+!   EXPLOSIVE: Apply volcanic emissions
+!   -----------------------------------
+    if (workspace%nVolcE > 0) then
+       if (associated(SO2EMVE)) SO2EMVE=0.0
+       allocate(iPoint(workspace%nVolcE), jPoint(workspace%nVolcE),  __STAT__)
+       call MAPL_GetHorzIJIndex(workspace%nVolcE, iPoint, jPoint, &
+                                grid = grid,               &
+                                lon  = workspace%vLon/real(MAPL_RADIANS_TO_DEGREES), &
+                                lat  = workspace%vLat/real(MAPL_RADIANS_TO_DEGREES), &
+                                rc   = status)
+           if ( status /= 0 ) then
+              if (mapl_am_i_root()) print*, trim(Iam), ' - cannot get indices for point emissions'
+              VERIFY_(status)
+           end if
+
+       call SUvolcanicEmissions (workspace%nVolcE, workspace%vStartE, workspace%vEndE, workspace%vSO2E, workspace%vElevE, &
+                                 workspace%vCloudE, iPoint, jPoint, nhms, SO2EMVE, SO2, nSO2, SUEM, &
+                                 self%km, self%cdt, MAPL_GRAV, zle, delp, area, workspace%vLatE, workspace%vLonE, __RC__)
+       deallocate(iPoint, jPoint, __STAT__)
     end if
 
 !   Apply diurnal cycle if so desired
@@ -926,6 +978,7 @@ subroutine Run1 (GC, import, export, clock, RC)
                                        workspace%pStart, workspace%pEnd, zle, &
                                        area, iPoint, jPoint, nhms, emissions_point, __RC__)
 
+        deallocate(iPoint, jPoint, __STAT__)
         SO4 = SO4 + self%cdt * MAPL_GRAV / delp * emissions_point
      end if
 

ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_instance_SU.rc

@@ -8,7 +8,8 @@ aerosol_monochromatic_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opti
 nbins: 4
 
 # Volcanic pointwise sources
-volcano_srcfilen: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_CARN_v202106.degassing_only.rc
+volcano_srcfilen_explosive: /dev/null
+volcano_srcfilen_degassing: ExtData/chemistry/CARN/v202106/sfc/so2_volcanic_emissions_CARN_v202106.degassing_only.rc
 
 # Heights [m] of LTO, CDS and CRS aviation emissions layers
 aviation_vertical_layers: 0.0 100.0 9.0e3 10.0e3

Process_Library/GOCART2G_Process.F90

@@ -5673,7 +5673,7 @@ end subroutine DMSemission
 ! !IROUTINE: SUvolcanicEmissions
 
    subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint, &
-                                   jPoint, nhms, SO2EMVN, SO2EMVE, SO2, nSO2, SU_emis, km, cdt, grav,&
+                                   jPoint, nhms, SO2EMVol, SO2, nSO2, SU_emis, km, cdt, grav,&
                                    hghte, delp, area, vLat, vLon, rc)
 ! !USES:
    implicit NONE
@@ -5697,8 +5697,7 @@ subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint
    real, dimension(:), intent(in)     :: vLat  ! latitude specified in file [degree]
    real, dimension(:), intent(in)     :: vLon  ! longitude specified in file [degree]
 ! !INOUT PARAMETERS:
-  real, pointer, dimension(:,:), intent(inout) :: SO2EMVN ! non-explosive volcanic emissions [kg m-2 s-1]
-  real, pointer, dimension(:,:), intent(inout) :: SO2EMVE ! explosive volcanic emissions [kg m-2 s-1]
+  real, pointer, dimension(:,:), intent(inout) :: SO2EMVol ! volcanic emissions [kg m-2 s-1]
   real, pointer, dimension(:,:,:), intent(inout) :: SO2 ! SO2 [kg kg-1]
   real, pointer, dimension(:,:,:), intent(inout) :: SU_emis      ! SU emissions, kg/m2/s
   real, dimension(:), intent(inout) ::  vElev ! bottom elevation of emissions [m]
@@ -5722,7 +5721,6 @@ subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint
    real :: deltaSO2v
    real, dimension(:,:), allocatable :: z0
    real, allocatable, dimension(:,:) :: srcSO2volc
-   real, allocatable, dimension(:,:) :: srcSO2volce
 
 !EOP
 !-------------------------------------------------------------------------
@@ -5731,13 +5729,9 @@ subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint
    if (nVolc > 0) then
 
    allocate(srcSO2volc, mold=area)
-   allocate(srcSO2volce, mold=area)
    srcSO2volc = 0.
-   srcSO2volce = 0.
 
-   if (associated(SU_emis)) SU_emis = 0.0
-   if (associated(SO2EMVN)) SO2EMVN = 0.
-   if (associated(SO2EMVE)) SO2EMVE = 0.
+   if (associated(SO2EMVol)) SO2EMVol = 0.
 
    allocate(z0, mold=area)
    z0 = hghte(:,:,km)
@@ -5776,11 +5770,7 @@ subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint
 
 !        Diagnostic - sum of volcanos
 !        ----------------------------
-         if (hup .eq. hlow) then
-            srcSO2volc(i,j) = srcSO2volc(i,j) + so2volcano
-         else
-            srcSO2volce(i,j) = srcSO2volce(i,j) + so2volcano
-         endif
+         srcSO2volc(i,j) = srcSO2volc(i,j) + so2volcano
 
          dzvolc = hup-hlow
          do k = km, 1, -1
@@ -5827,9 +5817,8 @@ subroutine SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, vElev, vCloud, iPoint
    enddo     ! it
   end if ! nVolc > 0
 
-  if (associated(SO2EMVN)) SO2EMVN = SO2EMVN + srcSO2volc
-  if (associated(SO2EMVE)) SO2EMVE = SO2EMVE + srcSO2volce
-  if (associated(SU_emis)) SU_emis(:,:,nSO2) = SU_emis(:,:,nSO2) + srcSO2volc + srcSO2volce
+  if (associated(SO2EMVol)) SO2EMVol = SO2EMVol + srcSO2volc
+  if (associated(SU_emis)) SU_emis(:,:,nSO2) = SU_emis(:,:,nSO2) + srcSO2volc
 
   __RETURN__(__SUCCESS__)
   end subroutine SUvolcanicEmissions