Merge pull request NCAR#64 from SamuelTrahanNOAA/feature/v2-ca-pert-g…

…bbepx

SPPT and CA perturbations of SO2 emissions + 2 regtests

gsdchem/gsd_chem_plume_wrapper.F90

@@ -20,7 +20,19 @@ module gsd_chem_plume_wrapper
 
 !> \brief Brief description of the subroutine
 !!
-      subroutine gsd_chem_plume_wrapper_init()
+      subroutine gsd_chem_plume_wrapper_init(ca_global_emis,do_sppt_emis,im,emis_multiplier,errmsg,errflg)
+        implicit none
+        logical, intent(in) :: ca_global_emis, do_sppt_emis
+        real, intent(out) :: emis_multiplier(:)
+        character(len=*), intent(out) :: errmsg
+        integer, intent(out) :: errflg, im
+
+        errmsg=''
+        errflg=0
+
+        if(ca_global_emis .or. do_sppt_emis) then
+          emis_multiplier=1.0
+        endif
       end subroutine gsd_chem_plume_wrapper_init
 
 !> \brief Brief description of the subroutine
@@ -46,6 +58,7 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
                    ntrac,ntso2,ntpp25,ntbc1,ntoc1,ntpp10,                        &
                    gq0,qgrs,ebu,abem,                                            &
                    biomass_burn_opt_in,plumerise_flag_in,plumerisefire_frq_in,   &
+                   emis_multiplier, ca1, ca_global_emis, do_sppt_emis, sppt_wts, &
                    errmsg,errflg)
 
     implicit none
@@ -59,9 +72,13 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
     integer, parameter :: ims=1,jms=1,jme=1, kms=1
     integer, parameter :: its=1,jts=1,jte=1, kts=1
 
+    logical,        intent(in) :: ca_global_emis, do_sppt_emis
+    real, optional, intent(inout) :: emis_multiplier(:)
+    real, intent(in)    :: ca1(im)
     integer, dimension(im), intent(in) :: vegtype    
     real(kind_phys), dimension(im,    5), intent(in) :: fire_GBBEPx
     real(kind_phys), dimension(im,   13), intent(in) :: fire_MODIS
+    real(kind_phys), optional, intent(in) :: sppt_wts(:,:)
     real(kind_phys), dimension(im,kme), intent(in) :: ph3d, pr3d
     real(kind_phys), dimension(im,kte), intent(in) :: phl3d, prl3d, tk3d,        &
                 us3d, vs3d, spechum, w
@@ -97,8 +114,8 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
     real(kind_phys), parameter :: ugkg = 1.e-09_kind_phys !lzhang
 
 !>-- local variables
-    real(kind_phys) :: curr_secs
-    real(kind_phys) :: factor, factor2, factor3
+    real(kind_phys) :: curr_secs, ca1_scaled
+    real(kind_phys) :: factor, factor2, factor3, random_factor(ims:im)
     integer :: nbegin
     integer :: i, j, jp, k, kp, n
 
@@ -109,7 +126,7 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
     biomass_burn_opt  = biomass_burn_opt_in
     plumerise_flag    = plumerise_flag_in
     plumerisefire_frq = plumerisefire_frq_in
-
+    random_factor = 1.0
     curr_secs = ktau * dt
 
     ! -- set domain
@@ -171,6 +188,27 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
     ! -- add biomass burning emissions at every timestep
     if (biomass_burn_opt == BURN_OPT_ENABLE) then
       jp = jte
+
+      if(plumerise_flag == FIRE_OPT_GBBEPx) then
+       if (do_sppt_emis) then
+        do i = ims, im
+          emis_multiplier(i) = max(0.5,min(1.5,sppt_wts(i,kme/2)))
+          random_factor(i) = emis_multiplier(i)
+        enddo
+       elseif (ca_global_emis) then
+        do i = ims, im
+          ! ca1(i) is always precisely 0 or 2
+          if(ca1(i)<1.0) then
+            ca1_scaled=0.9
+          else
+            ca1_scaled=1.0/0.9
+          endif
+          emis_multiplier(i) = max(0.5,min(1.5,emis_multiplier(i)*0.95 + ca1_scaled*0.05))
+          random_factor(i) = emis_multiplier(i)
+        enddo
+       endif
+      endif
+
       factor3 = 0._kind_phys
       select case (plumerise_flag)
         case (FIRE_OPT_MODIS)
@@ -196,6 +234,9 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
             ! -- factor for pm emissions, factor2 for burn emissions
             factor  = dt*rri(i,k,j)/dz8w(i,k,j)
             factor2 = factor * factor3
+              if(plumerise_flag==FIRE_OPT_GBBEPx) then
+                factor2 = factor2 * random_factor(i)
+              endif
             chem(i,k,j,p_oc1) = chem(i,k,j,p_oc1) + factor  * ebu_in(i,j,p_ebu_in_oc  )
             chem(i,k,j,p_bc1) = chem(i,k,j,p_bc1) + factor  * ebu_in(i,j,p_ebu_in_bc  )
             chem(i,k,j,p_p25) = chem(i,k,j,p_p25) + factor  * ebu_in(i,j,p_ebu_in_pm25)
@@ -212,6 +253,9 @@ subroutine gsd_chem_plume_wrapper_run(im, kte, kme, ktau, dt,                &
               ! -- factor for pm emissions, factor2 for burn emissions
               factor  = dt*rri(i,k,j)/dz8w(i,k,j)
               factor2 = factor * factor3
+              if(plumerise_flag==FIRE_OPT_GBBEPx) then
+                factor2 = factor2 * random_factor(i)
+              endif
               chem(i,k,j,p_oc1) = chem(i,k,j,p_oc1) + factor  * ebu(i,k,j,p_ebu_oc  )
               chem(i,k,j,p_bc1) = chem(i,k,j,p_bc1) + factor  * ebu(i,k,j,p_ebu_bc  )
               chem(i,k,j,p_p25) = chem(i,k,j,p_p25) + factor  * ebu(i,k,j,p_ebu_pm25)

gsdchem/gsd_chem_plume_wrapper.meta

@@ -7,6 +7,56 @@
 [ccpp-arg-table]
   name = gsd_chem_plume_wrapper_init
   type = scheme
+[ca_global_emis]
+  standard_name = flag_for_tracer_emissions_global_cellular_automata
+  long_name = switch for global ca applied to chemical tracer emissions
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[do_sppt_emis]
+  standard_name = flag_for_stochastic_emissions_perturbations
+  long_name = flag for stochastic emissions perturbations
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[im]
+  standard_name = horizontal_dimension
+  long_name = horizontal dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+  optional = F
+[emis_multiplier]
+  standard_name = gsd_chem_ca_global_emis_multiplier
+  long_name = fraction of emissions to generate based on cellular automata
+  units = frac
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
+  optional = F
+[errmsg]
+  standard_name = ccpp_error_message
+  long_name = error message for error handling in CCPP
+  units = none
+  dimensions = ()
+  type = character
+  kind = len=*
+  intent = out
+  optional = F
+[errflg]
+  standard_name = ccpp_error_flag
+  long_name = error flag for error handling in CCPP
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = out
+  optional = F
 
 ########################################################################
 [ccpp-arg-table]
@@ -273,6 +323,49 @@
   type = integer
   intent = in
   optional = F
+[emis_multiplier]
+  standard_name = gsd_chem_ca_global_emis_multiplier
+  long_name = fraction of emissions to generate based on cellular automata
+  units = frac
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
+[ca1]
+  standard_name = cellular_automata_global_pattern
+  long_name = cellular automata global pattern
+  units = flag
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[ca_global_emis]
+  standard_name = flag_for_tracer_emissions_global_cellular_automata
+  long_name = switch for global ca applied to chemical tracer emissions
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[do_sppt_emis]
+  standard_name = flag_for_stochastic_emissions_perturbations
+  long_name = flag for stochastic emissions perturbations
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[sppt_wts]
+  standard_name = weights_for_stochastic_sppt_perturbation
+  long_name = weights for stochastic sppt perturbation
+  units = none
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP