[FEATURE REQUEST] Add emissions to State_Chm%Species and subtract drydep from State_Chm%Species before PBL mixing

[yantosca]

In order to further clean up the PBL mixing schemes, we will need to do the following:

1. Call HEMCO to update Emissions fluxes
2. Add emissions fluxes into State_Chm%Species (all levels)
3. Call DO_DRYDEP to update dry-deposition frequencies.
4. Remove the dry-deposition contribution from State_Chm%Species (surface level only)
5. Call either VDIFF (non-local PBL mixing) or TURBDAY (full PBL mixing)

Adding the emissions and subtracting the drydep contributions separately will also allow us to archive budget terms for just emissions and just drydep.

Currently, the code applies a surface flux (Emissions - Drydep) within the non-local PBL mixing scheme. This is problematic because calls to HEMCO code was originally embedded within the non-local PBL scheme (but this has now been removed). Also, the current setup can only let us compute a budget term for Emissions+Drydep instead of Emissions and Drydep separately.

[lizziel]

@yantosca, can you give follow-up on this work here? Is it all done?

[yantosca]

So far I was only able to move the HEMCO-specific code out of vdiff_mod.F90 and into hco_utilities_mod.F90 in routine Compute_Sflx_for_Vdiff. This stores the surface flux (formerly sflx in VDIFFDR) into State_Chm%SurfaceFlux. This was the fastest thing to do. That was what I thought the scope of the problem was. Perhaps I misunderstood it?

[lizziel]

@jimmielin, could you expand on the remainder of work that is needed and by when? Then we can hone on in what milestone to include this in and make sure we have a full list of what needs to be done. Perhaps some of it can go into 13.0, or all of it can go into 13.1.

[yantosca]

Also -- we are prioritizing 13.0.0 items for now, since we need to get that out later this summer (or early fall). I think that to implemetn the full removal of Emissions & drydep from mixing might result in numerical differences. For 13.0.0 we need to include only updates that won't change the results that much (i.e. no science or structural updates that drastically change output).

[jimmielin]

Hi Bob, may I ask how the 3-D emissions (above surface) are being handled? I could not locate where the emissions within the PBL but above the surface are handled, as State_Chm%SurfaceFlux is a (I, J, N) "2-D" array.

It appears that the other place where emissions are applied, in Mixing_Mod.F90:Do_Tend, only handles emissions above the PBL when non-local PBL mixing is used:

  SUBROUTINE Do_Mixing( Input_Opt,  State_Chm, State_Diag, &
                        State_Grid, State_Met, RC )
...
    !-----------------------------------------------------------------------
    ! Do non-local PBL mixing. This will apply the species tendencies
    ! (emission fluxes and dry deposition rates) below the PBL.
    ! This is done for all species with defined emissions / dry
    ! deposition rates, including dust.
    !
    ! Set OnlyAbovePBL flag (used below by DO_TEND) to indicate that
    ! fluxes within the PBL have already been applied.
    ! ----------------------------------------------------------------------
    IF ( Input_Opt%LTURB .AND. Input_Opt%LNLPBL ) THEN

       ! Non-local mixing
       CALL Do_Vdiff( Input_Opt,  State_Chm, State_Diag,                     &
                      State_Grid, State_Met, RC                             )
...
       ! Fluxes in PBL have been applied (non-local PBL mixing)
       OnlyAbovePBL = .TRUE.

    ELSE

       ! Fluxes in PBL have not been applied (full PBL mixing)
       OnlyAbovePBL = .FALSE.

    ENDIF
...
    !-----------------------------------------------------------------------
    ! Apply tendencies. This will apply dry deposition rates and
    ! emission fluxes below the PBL if it has not yet been done
    ! via the non-local PBL mixing. It also adds the emissions above
    ! the PBL to the species array. Emissions of some species may be
    ! capped at the tropopause to avoid build-up in stratosphere.
    !-----------------------------------------------------------------------

    ! Apply tendencies
    CALL DO_TEND( Input_Opt, State_Chm,  State_Diag, &
                  State_Grid, State_Met, OnlyAbovePBL, RC )

Since here OnlyAbovePBL = .true. I am a little puzzled as to where the 3-D emissions for levels >= 2 and below the PBL level are applied. Sorry I couldn't find them, could you please clarify? Thanks!

[yantosca]

I agree this is very confusing...this goes back to Christoph's original implementation of DO_TEND.

In the old vdiff_mod.F90 (and now in the new Compute_Sflx_for_Vdiff) the emissions within the PBL are summed together and stored in State_Chm%SurfaceFlux. This is done starting at line 1712 below:

geos-chem/GeosCore/hco_utilities_gc_mod.F90

Lines 1669 to 1721 in 17bc5f8

	DO J = 1, State_Grid%NY
	DO I = 1, State_Grid%NX
	! PBL top level [integral model levels]
	topMix = MAX( 1, FLOOR( State_Met%PBL_TOP_L(I,J) ) )
	! Loop over advected species
	DO NA = 1, State_Chm%nAdvect
	! Get the modelId
	N = State_Chm%Map_Advect(NA)
	! Point to the corresponding entry in the species database
	ThisSpc => State_Chm%SpcData(N)%Info
	!------------------------------------------------------------------
	! Add total emissions in the PBL to the EFLX array
	! which tracks emission fluxes. Units are [kg/m2/s].
	!------------------------------------------------------------------
	IF ( Input_Opt%ITS_A_CH4_SIM ) THEN
	! CH4 emissions become stored in CH4_EMIS in global_ch4_mod.F90.
	! We use CH4_EMIS here instead of the HEMCO internal emissions
	! only to make sure that total CH4 emissions are properly defined
	! in a multi-tracer CH4 simulation. For a single-tracer simulation
	! and/or all other source types, we could use the HEMCO internal
	! values set above and would not need the code below.
	! Units are already in kg/m2/s. (ckeller, 10/21/2014)
	!
	!%%% NOTE: MAYBE THIS CAN BE REMOVED SOON (bmy, 5/18/19)%%%
	eflx(I,J,NA) = CH4_EMIS(I,J,NA)
	ELSE IF ( Input_Opt%ITS_A_MERCURY_SIM ) THEN
	! HG emissions become stored in HG_EMIS in mercury_mod.F90.
	! This is a workaround to ensure backwards compatibility.
	! Units are already in kg/m2/s. (ckeller, 10/21/2014)
	!
	!%%% NOTE: MAYBE THIS CAN BE REMOVED SOON (bmy, 5/18/19)%%%
	eflx(I,J,NA) = HG_EMIS(I,J,NA)
	ELSE
	! Compute emissions for all other simulation
	tmpFlx = 0.0_fp
	DO L = 1, topMix
	CALL GetHcoValEmis( NA, I, J, L, found, emis )
	IF ( .NOT. found ) EXIT
	tmpFlx = tmpFlx + emis
	ENDDO
	eflx(I,J,NA) = eflx(I,J,NA) + tmpFlx
	ENDIF

All of the emissions go into the surface layer (hence sflx is the variable name used into vdiff_mod.F90). Then the VDIFF routine partitions that into the PBL vertically.

[jimmielin]

Thanks Bob! I agree this is very confusing. We should probably confirm that distributing emissions in the PBL like this is no longer necessary. I spoke with Daniel a few months back and he mentioned this was originally added because of limited vertical levels many many versions before, resulting in "lumping" of emissions in the surface layer, and thus distributing in the PBL was necessary.

Once the scientific implications are sorted out, the to-do items would be:

• Get rid of the code to perform emission fluxes from PBL mixing;
• Store emission fluxes in a (name TBD) State_Chm%Flux 3-D array (I, J, K, N);
• Apply all tendencies in DO_TEND over the entire vertical, regardless of PBL mixing method.

This would allow input of emissions and dry deposition frequencies from a module external to the GEOS-Chem core, through this new flux array in State_Chm. This will be most relevant for CESM2-GC.