Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bug fixes for the Hg simulation with KPP #1336

Merged
merged 7 commits into from Aug 8, 2022
Merged

Bug fixes for the Hg simulation with KPP #1336

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
b438531
Bug fix: replace VAR with C in Prod/Loss diag section of fullchem_mod…
yantosca Jun 21, 2022
3c676ec
Update species_database.yml for Hg simulation following Shah et al 2021
yantosca Jun 23, 2022
4c8f6cf
Update species_database_hg.yml to add Is_Hg2 and Is_HgP tags
yantosca Jun 23, 2022
2c4e940
Re-add MercuryEmis collection to HISTORY.rc.Hg template file
yantosca Jun 23, 2022
2792f85
Restore ocean Hg and snow Hg quantities to HISTORY.rc.Hg template file
yantosca Jun 23, 2022
e50fc2f
Remove all references to N_HG_CATS and related code from Hg routines
yantosca Jun 23, 2022
10003b8
Add missing call to SrcHg0 in routine EmissMercury
yantosca Jun 29, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
14 changes: 4 additions & 10 deletions GeosCore/convection_mod.F90
View file
Open in desktop
Expand Up @@ -1290,12 +1290,9 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt, &
! Wet scavenged Hg(II) in [kg]
WET_Hg2 = ( T0_SUM * AREA_M2 )

! Category # for this Hg2 species
Hg_Cat = SpcInfo%Hg_Cat

! Pass to "ocean_mercury_mod.f"
CALL ADD_Hg2_WD ( I, J, Hg_Cat, WET_Hg2 )
CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, WET_Hg2, &
CALL ADD_Hg2_WD ( I, J, WET_Hg2 )
CALL ADD_Hg2_SNOWPACK( I, J, WET_Hg2, &
State_Met, State_Chm, State_Diag )
ENDIF

Expand All @@ -1307,12 +1304,9 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt, &
! Wet scavenged Hg(P) in [kg]
WET_HgP = ( T0_SUM * AREA_M2 )

! Category # for this Hg2 species
Hg_Cat = SpcInfo%Hg_Cat

! Pass to "ocean_mercury_mod.f"
CALL ADD_HgP_WD ( I, J, Hg_Cat, WET_HgP )
CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, WET_HgP, &
CALL ADD_HgP_WD ( I, J, WET_HgP )
CALL ADD_Hg2_SNOWPACK( I, J, WET_HgP, &
State_Met, State_Chm, State_Diag )
ENDIF
ENDIF
Expand Down
159 changes: 37 additions & 122 deletions GeosCore/depo_mercury_mod.F90
View file
Open in desktop
Expand Up @@ -42,10 +42,10 @@ MODULE DEPO_MERCURY_MOD
PUBLIC :: HG2mth_wd, HG0mth_dd, HG2mth_dd
PUBLIC :: LHGSNOW

REAL(fp), ALLOCATABLE :: DD_Hg2(:,:,:)
REAL(fp), ALLOCATABLE :: DD_HgP(:,:,:)
REAL(fp), ALLOCATABLE :: WD_Hg2(:,:,:)
REAL(fp), ALLOCATABLE :: WD_HgP(:,:,:)
REAL(fp), ALLOCATABLE :: DD_Hg2(:,:)
REAL(fp), ALLOCATABLE :: DD_HgP(:,:)
REAL(fp), ALLOCATABLE :: WD_Hg2(:,:)
REAL(fp), ALLOCATABLE :: WD_HgP(:,:)
REAL(fp), ALLOCATABLE :: HG0mth_dd(:,:)
REAL(fp), ALLOCATABLE :: HG2mth_dd(:,:)
REAL(fp), ALLOCATABLE :: HG2mth_wd(:,:)
Expand All @@ -68,18 +68,8 @@ MODULE DEPO_MERCURY_MOD
! !PRIVATE TYPES:
!
! Scalars for Hg indexing
INTEGER :: N_Hg_CATS
INTEGER :: id_Hg0
INTEGER :: id_Hg2

! Pointers for Hg Indexing
! NOTE: Because these are SAVEd pointers (by virtue of being
! module variables) we can set these to NULL here. (bmy, 4/29/16)
INTEGER, POINTER :: Hg0_Id_List (:) => NULL()
INTEGER, POINTER :: Hg2_Id_List (:) => NULL()
INTEGER, POINTER :: HgP_Id_List (:) => NULL()
INTEGER, POINTER :: Hg0_CAT(:) => NULL()
INTEGER, POINTER :: Hg2_CAT(:) => NULL()
INTEGER :: id_Hg0
INTEGER :: id_Hg2

CONTAINS
!EOC
Expand All @@ -96,12 +86,11 @@ MODULE DEPO_MERCURY_MOD
!\\
! !INTERFACE:
!
SUBROUTINE ADD_Hg2_DD( I, J, HG_CAT, DRY_Hg2 )
SUBROUTINE ADD_Hg2_DD( I, J, DRY_Hg2 )
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J ! Grid box lon & lat indices
INTEGER, INTENT(IN) :: Hg_Cat ! Hg category number
REAL(fp), INTENT(IN) :: DRY_Hg2 ! Hg(II) dry deposited out of the
! atmosphere [kg]
!
Expand All @@ -113,9 +102,7 @@ SUBROUTINE ADD_Hg2_DD( I, J, HG_CAT, DRY_Hg2 )
!BOC

! Store dry deposited Hg(II) into DD_Hg2 array
IF ( Hg_Cat > 0 ) THEN
DD_Hg2(I,J,Hg_Cat) = DD_Hg2(I,J,Hg_Cat) + DRY_Hg2
ENDIF
DD_Hg2(I,J) = DD_Hg2(I,J) + DRY_Hg2

END SUBROUTINE ADD_Hg2_DD
!EOC
Expand All @@ -132,12 +119,11 @@ END SUBROUTINE ADD_Hg2_DD
!\\
! !INTERFACE:
!
SUBROUTINE ADD_Hg2_WD( I, J, Hg_Cat, WET_Hg2 )
SUBROUTINE ADD_Hg2_WD( I, J, WET_Hg2 )
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J ! Grid box lon & lat indices
INTEGER, INTENT(IN) :: Hg_Cat ! Hg category number
REAL(fp), INTENT(IN) :: WET_Hg2 ! Hg(II) scavenged out of the
! atmosphere [kg]
!
Expand All @@ -149,9 +135,7 @@ SUBROUTINE ADD_Hg2_WD( I, J, Hg_Cat, WET_Hg2 )
!BOC

! Store wet deposited Hg(II) into WD_Hg2 array
IF ( Hg_Cat > 0 ) THEN
WD_Hg2(I,J,Hg_Cat) = WD_Hg2(I,J,Hg_Cat) + WET_Hg2
ENDIF
WD_Hg2(I,J) = WD_Hg2(I,J) + WET_Hg2

END SUBROUTINE ADD_Hg2_WD
!EOC
Expand All @@ -168,12 +152,11 @@ END SUBROUTINE ADD_Hg2_WD
!\\
! !INTERFACE:
!
SUBROUTINE ADD_HgP_DD( I, J, Hg_Cat, DRY_HgP )
SUBROUTINE ADD_HgP_DD( I, J, DRY_HgP )
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J ! Grid box lon & lat indices
INTEGER, INTENT(IN) :: Hg_Cat ! Hg category number
REAL(fp), INTENT(IN) :: DRY_HgP ! HgP dry deposited out of the
! atmosphere [kg]
!
Expand All @@ -185,9 +168,7 @@ SUBROUTINE ADD_HgP_DD( I, J, Hg_Cat, DRY_HgP )
!BOC

! Store dry deposited Hg(II) into DD_Hg2 array
IF ( Hg_Cat > 0 ) THEN
DD_HgP(I,J,Hg_Cat) = DD_HgP(I,J,Hg_Cat) + DRY_HgP
ENDIF
DD_HgP(I,J) = DD_HgP(I,J) + DRY_HgP

END SUBROUTINE ADD_HgP_DD
!EOC
Expand All @@ -204,12 +185,11 @@ END SUBROUTINE ADD_HgP_DD
!\\
! !INTERFACE:
!
SUBROUTINE ADD_HgP_WD( I, J, Hg_Cat, WET_HgP )
SUBROUTINE ADD_HgP_WD( I, J, WET_HgP )
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J ! Grid box lon & lat indices
INTEGER, INTENT(IN) :: Hg_Cat ! Hg category number
REAL(fp), INTENT(IN) :: WET_HgP ! HgP scavenged out of the
! atmosphere [kg]
!
Expand All @@ -221,9 +201,7 @@ SUBROUTINE ADD_HgP_WD( I, J, Hg_Cat, WET_HgP )
!BOC

! Store wet deposited HgP into WD_HgP array
IF ( Hg_Cat > 0 ) THEN
WD_HgP(I,J,Hg_Cat) = WD_HgP(I,J,Hg_Cat) + WET_HgP
ENDIF
WD_HgP(I,J) = WD_HgP(I,J) + WET_HgP

END SUBROUTINE ADD_HgP_WD
!EOC
Expand All @@ -239,8 +217,7 @@ END SUBROUTINE ADD_HgP_WD
!\\
! !INTERFACE:
!
SUBROUTINE ADD_HG2_SNOWPACK( I, J, Hg_Cat, DEP_Hg2, &
State_Met, State_Chm, State_Diag )
SUBROUTINE ADD_HG2_SNOWPACK( I, J, DEP_Hg2, State_Met, State_Chm, State_Diag )
!
! !USES:
!
Expand All @@ -252,7 +229,6 @@ SUBROUTINE ADD_HG2_SNOWPACK( I, J, Hg_Cat, DEP_Hg2, &
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J ! Grid box lon & lat indices
INTEGER, INTENT(IN) :: Hg_Cat ! Hg category number
REAL(fp), INTENT(IN) :: Dep_Hg2 ! Hg2 (or HgP) deposited
TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object
!
Expand All @@ -270,17 +246,16 @@ SUBROUTINE ADD_HG2_SNOWPACK( I, J, Hg_Cat, DEP_Hg2, &
!
! !LOCAL VARIABLES:
!
INTEGER :: NN
LOGICAL :: IS_SNOW_OR_ICE
REAL(fp) :: FRAC_SNOW_OR_ICE
REAL(fp) :: FRAC_O
REAL(fp) :: DT

! Pointers
REAL(fp), POINTER :: SNOW_HG_OC(:,:,:)
REAL(fp), POINTER :: SNOW_HG_LN(:,:,:)
REAL(fp), POINTER :: SNOW_HG_STORED_OC(:,:,:)
REAL(fp), POINTER :: SNOW_HG_STORED_LN(:,:,:)
REAL(fp), POINTER :: SNOW_HG_OC(:,:)
REAL(fp), POINTER :: SNOW_HG_LN(:,:)
REAL(fp), POINTER :: SNOW_HG_STORED_OC(:,:)
REAL(fp), POINTER :: SNOW_HG_STORED_LN(:,:)

!=================================================================
! ADD_HG2_SNOWPACK begins here!
Expand All @@ -292,9 +267,6 @@ SUBROUTINE ADD_HG2_SNOWPACK( I, J, Hg_Cat, DEP_Hg2, &
SNOW_HG_STORED_OC => State_Chm%SnowHgOceanStored
SNOW_HG_STORED_LN => State_Chm%SnowHgLandStored

! Store the Hg category number in NN
NN = Hg_Cat

! Return if snowpack model is disabled
IF ( .NOT. LHGSNOW ) RETURN

Expand All @@ -316,27 +288,27 @@ SUBROUTINE ADD_HG2_SNOWPACK( I, J, Hg_Cat, DEP_Hg2, &
! Add 60% of deposition to surface (i.e. assume 40% accumulates
! in surface) multiplied by the fraction of the box that is
! snow or ice (i.e. 1 for all met fields but MERRA and GEOS-5.7.x)
SNOW_HG_OC(I,J,NN) = SNOW_HG_OC(I,J,NN) + &
FRAC_O * FRAC_SNOW_OR_ICE * &
MAX( 0.6e+0_fp*DEP_HG2, 0e+0_fp )
SNOW_HG_OC(I,J) = SNOW_HG_OC(I,J) + &
FRAC_O * FRAC_SNOW_OR_ICE * &
MAX( 0.6e+0_fp*DEP_HG2, 0e+0_fp )

! Add remaining deposited Hg to reservoir for delivery to ocean
! when snow melts later (jaf, 6/17/11)
! This is Hg in snowpack over ocean that CANNOT be
! re-emitted to atmosphere
SNOW_HG_STORED_OC(I,J,NN) = SNOW_HG_STORED_OC(I,J,NN) + &
FRAC_O * FRAC_SNOW_OR_ICE * &
MAX( 0.4e+0_fp*DEP_HG2, 0e+0_fp )
SNOW_HG_STORED_OC(I,J) = SNOW_HG_STORED_OC(I,J) + &
FRAC_O * FRAC_SNOW_OR_ICE * &
MAX( 0.4e+0_fp*DEP_HG2, 0e+0_fp )

! This is Hg in snowpack over land that can potentially be
! re-emitted to atmosphere
SNOW_HG_LN(I,J,NN) = SNOW_HG_LN(I,J,NN) + &
(1e+0_fp - FRAC_O) * FRAC_SNOW_OR_ICE * &
MAX( 0.6e+0_fp*DEP_HG2, 0e+0_fp )
SNOW_HG_LN(I,J) = SNOW_HG_LN(I,J) + &
(1e+0_fp - FRAC_O) * FRAC_SNOW_OR_ICE * &
MAX( 0.6e+0_fp*DEP_HG2, 0e+0_fp )

! This is Hg in snowpack over land that CANNOT be
! re-emitted to atmosphere
SNOW_HG_STORED_LN(I,J,NN) = SNOW_HG_STORED_LN(I,J,NN) + &
SNOW_HG_STORED_LN(I,J) = SNOW_HG_STORED_LN(I,J) + &
( 1e+0_fp - FRAC_O ) * &
FRAC_SNOW_OR_ICE &
* MAX( 0.4e+0_fp*DEP_HG2, 0e+0_fp )
Expand Down Expand Up @@ -926,79 +898,30 @@ SUBROUTINE INIT_DEPO_MERCURY( Input_Opt, State_Chm, State_Grid, RC )
! Exit if this is a dry-run
IF ( Input_Opt%DryRun ) RETURN

! GTMM restart file name
GTMM_RST_FILE = Input_Opt%GTMM_RST_FILE

!================================================================
! Initialize local Hg indexing variables
!================================================================

! Now get # of tagHg categories from State_Chm
N_Hg_CATS = State_Chm%N_Hg_CATS

! Hg species index corresponding to a given Hg category number
Hg0_Id_List => State_Chm%Hg0_Id_List
Hg2_Id_List => State_Chm%Hg2_Id_List
HgP_Id_List => State_Chm%HgP_Id_List

! Category numbers for each Hg0 tracer
ALLOCATE( Hg0_Cat( State_Chm%nSpecies ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'Hg0_Cat' )
Hg0_Cat = 0

! Category numbers for each Hg2 tracer
ALLOCATE( Hg2_Cat( State_Chm%nSpecies ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'Hg0_Cat' )
Hg2_Cat = 0

! Loop ovver all Hg species
DO N = 1, State_Chm%nSpecies

! Point to the Species Database entry for tracer # N
SpcInfo => State_Chm%SpcData(N)%Info

! Define local id_Hg0 and id_Hg2 flags
SELECT CASE( TRIM( SpcInfo%Name ) )
CASE( 'Hg0' )
id_Hg0 = SpcInfo%ModelId
CASE( 'Hg2' )
id_Hg2 = SpcInfo%ModelId
CASE DEFAULT
! Do nothing
END SELECT

! Store the Hg0 and Hg2 category numbers for each tracer
! for future lookup in READ_GTMM_RESTART (bmy, 4/26/16)
IF ( SpcInfo%Is_Hg0 ) THEN
Hg0_Cat(N) = SpcInfo%Hg_Cat
ELSE iF ( SpcInfo%Is_Hg2 ) THEN
Hg2_Cat(N) = SpcInfo%Hg_Cat
ENDIF

! Free pointer
SpcInfo => NULL()
ENDDO

!================================================================
! Allocate arrays
!================================================================
ALLOCATE( DD_Hg2( State_Grid%NX, State_Grid%NY, N_Hg_CATS ), STAT=RC )
ALLOCATE( DD_Hg2( State_Grid%NX, State_Grid%NY ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'DD_Hg2' )
DD_Hg2 = 0e+0_fp

ALLOCATE( DD_HgP( State_Grid%NX, State_Grid%NY, N_Hg_CATS ), STAT=RC )
ALLOCATE( DD_HgP( State_Grid%NX, State_Grid%NY ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'DD_HgP' )
DD_HgP = 0e+0_fp

ALLOCATE( WD_Hg2( State_Grid%NX, State_Grid%NY, N_Hg_CATS ), STAT=RC )
ALLOCATE( WD_Hg2( State_Grid%NX, State_Grid%NY ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'WD_Hg2' )
WD_Hg2 = 0e+0_fp

ALLOCATE( WD_HgP( State_Grid%NX, State_Grid%NY, N_Hg_CATS ), STAT=RC )
ALLOCATE( WD_HgP( State_Grid%NX, State_Grid%NY ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'WD_HgP' )
WD_HgP = 0e+0_fp

IF ( Input_Opt%LGTMM ) THEN

! GTMM restart file name
GTMM_RST_FILE = Input_Opt%GTMM_RST_FILE

ALLOCATE( HG0mth_dd( State_Grid%NX, State_Grid%NY ), STAT=RC )
IF ( RC /= 0 ) CALL ALLOC_ERR( 'HG0mth_dd' )
HG0mth_dd = 0e+0_fp
Expand Down Expand Up @@ -1043,14 +966,6 @@ SUBROUTINE CLEANUP_DEPO_MERCURY
IF ( ALLOCATED( HG2mth_dd ) ) DEALLOCATE( HG2mth_dd )
IF ( ALLOCATED( HG2mth_wd ) ) DEALLOCATE( HG2mth_wd )

! Free pointers
Hg0_Id_List => NULL()
Hg2_Id_List => NULL()
HgP_Id_List => NULL()

IF ( ASSOCIATED( Hg0_Cat ) ) DEALLOCATE( Hg0_CAT )
IF ( ASSOCIATED( Hg2_Cat ) ) DEALLOCATE( Hg2_CAT )

END SUBROUTINE CLEANUP_DEPO_MERCURY
!EOC
END MODULE DEPO_MERCURY_MOD
14 changes: 14 additions & 0 deletions GeosCore/emissions_mod.F90
View file
Open in desktop
Expand Up @@ -149,6 +149,7 @@ SUBROUTINE Emissions_Run( Input_Opt, State_Chm, State_Diag, &
USE GLOBAL_CH4_MOD, ONLY : EMISSCH4
USE HCO_Interface_GC_Mod, ONLY : HCOI_GC_RUN
USE Input_Opt_Mod, ONLY : OptInput
USE Mercury_Mod, ONLY : EmissMercury
USE Precision_Mod
USE State_Chm_Mod, ONLY : ChmState
USE State_Diag_Mod, ONLY : DgnState
Expand Down Expand Up @@ -286,6 +287,19 @@ SUBROUTINE Emissions_Run( Input_Opt, State_Chm, State_Diag, &
ENDIF
ENDIF

! For mercury, use old emissions code for now
IF ( Input_Opt%ITS_A_MERCURY_SIM ) THEN
CALL EmissMercury( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, RC )

! Trap potential errors
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Error encountered in "EmissMercury"!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
ENDIF

! Prescribe some concentrations if needed
IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN
! Set other (non-UCX) fixed VMRs
Expand Down