Merge branch 'tangq/atm/chemUCI_surf_emission_fix' into NGD_v3atm (PR #3

)

Fix chemUCI high NOx

This PR corrects the surface emissions for chemUCI and fixes the high
NOx issue. The NO surface emissions could be added in upper troposphere
and lower stratosphere. With this fix, it will be limited to the
bottom layers.

The test has be done successfully and documented at this page.
The results show that the high NOx is gone.

[non-BFB]

components/eam/src/chemistry/mozart/mo_gas_phase_chemdr.F90

@@ -50,8 +50,8 @@ module mo_gas_phase_chemdr
   character(len=fieldname_len),dimension(extcnt)        :: extfrc_name
   logical :: convproc_do_aer 
 
-  real(r8), parameter :: low_limit = 0.1_r8 ! lower limit factor for dry deposition
-  real(r8), parameter :: up_limit  = 2.0_r8 ! upper limit factor for surface emission
+  real(r8), parameter :: low_limit = 0.1_r8  ! lower limit factor for dry deposition
+  integer,  parameter :: srf_emit_nlayer = 3 ! number of layers adding surface emission
 
 contains
 
@@ -494,7 +494,7 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, &
     real(r8) :: vmr0(ncol,pver,gas_pcnst)
     real(r8) :: vmr_old(ncol,pver,gas_pcnst)
     real(r8) :: vmr_old2(ncol,pver,gas_pcnst)
-    real(r8) :: tmp
+    real(r8) :: tmp(ncol)
 
     ! for explicit chem prodcution and loss
     real(r8)     ::  chem_prod(ncol,pver,gas_pcnst)         ! xported species (vmr/delt)
@@ -1545,37 +1545,21 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, &
     !-----------------------------------------------------------------------      
     !         ... Surface emissions
     !----------------------------------------------------------------------- 
-    ! if chemUCI is used, apply surface emissions here
+    ! if chemUCI is used, apply surface emissions to the bottom (srf_emit_nlayer) layers here
     if (uci1_ndx > 0) then 
        do m = 1,pcnst
           n = map2chm( m )
           if ( n > 0 ) then
             if ( any( cflx(:ncol,m) /= 0._r8 ) ) then
               if ( .not. any( aer_species == n ) .and. adv_mass(n) /= 0._r8 ) then
-                do k = pver, 1, -1 ! loop from bottom to top
+                do k = pver, pver-srf_emit_nlayer+1, -1 ! loop from bottom to top
                   ! kg/m2, tracer mass
                   wrk(:ncol,k) = adv_mass(n)*vmr(:ncol,k,n)/mbar(:ncol,k) &
                                     *pdeldry(:ncol,k)*rga
-                  j = 0 ! number of columns will double concentration after adding surf. emission
-                  do i = 1,ncol
-                    if ( cflx(i,m) /= 0._r8 ) then
-                      if ( wrk(i,k)*(up_limit-1._r8) >= cflx(i,m)*delt ) then
-                        tmp = wrk(i,k) + cflx(i,m)*delt
-                        cflx(i,m) = 0._r8
-                      else
-                        tmp = wrk(i,k)*up_limit
-                        cflx(i,m) = cflx(i,m) - wrk(i,k)*(up_limit-1._r8)*delt_inverse
-                        j = j + 1
-                      endif
-                      vmr(i,k,n) = tmp*mbar(i,k)/adv_mass(n)/pdeldry(i,k)/rga
-                    endif
-                  end do
-
-                  if ( j == 0 ) then
-                    exit
-                  endif
-
+                  tmp(:ncol) = wrk(:ncol,k) + cflx(:ncol,m)*delt/dble(srf_emit_nlayer)
+                  vmr(:ncol,k,n) = tmp(:ncol)*mbar(:ncol,k)/adv_mass(n)/pdeldry(:ncol,k)/rga
                 end do
+                cflx(:ncol,m) = 0._r8
               endif
             endif
           endif
@@ -1674,14 +1658,14 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, &
                   do i = 1,ncol
                     if ( sflx2(i,n) /= 0._r8 ) then
                       if ( wrk(i,k)*(1._r8-low_limit) >= sflx2(i,n)*delt ) then
-                        tmp = wrk(i,k) - sflx2(i,n)*delt
+                        tmp(1) = wrk(i,k) - sflx2(i,n)*delt
                         sflx2(i,n) = 0._r8
                       else
-                        tmp = wrk(i,k)*low_limit
+                        tmp(1) = wrk(i,k)*low_limit
                         sflx2(i,n) = sflx2(i,n) - wrk(i,k)*(1._r8-low_limit)*delt_inverse
                         j = j + 1
                       endif
-                      vmr(i,k,n) = tmp*mbar(i,k)/adv_mass(n)/pdeldry(i,k)/rga
+                      vmr(i,k,n) = tmp(1)*mbar(i,k)/adv_mass(n)/pdeldry(i,k)/rga
                     endif
                   end do