diff --git a/GeosCore/convection_mod.F90 b/GeosCore/convection_mod.F90 index a5b77a2f4..31877a5bc 100644 --- a/GeosCore/convection_mod.F90 +++ b/GeosCore/convection_mod.F90 @@ -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 @@ -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 diff --git a/GeosCore/depo_mercury_mod.F90 b/GeosCore/depo_mercury_mod.F90 index b764a26d6..3303e9ac9 100644 --- a/GeosCore/depo_mercury_mod.F90 +++ b/GeosCore/depo_mercury_mod.F90 @@ -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(:,:) @@ -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 @@ -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] ! @@ -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 @@ -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] ! @@ -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 @@ -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] ! @@ -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 @@ -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] ! @@ -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 @@ -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: ! @@ -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 ! @@ -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! @@ -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 @@ -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 ) @@ -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 @@ -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 diff --git a/GeosCore/emissions_mod.F90 b/GeosCore/emissions_mod.F90 index 9b6bb9bab..d71f54596 100644 --- a/GeosCore/emissions_mod.F90 +++ b/GeosCore/emissions_mod.F90 @@ -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 @@ -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 diff --git a/GeosCore/fullchem_mod.F90 b/GeosCore/fullchem_mod.F90 index 295a8c3fc..9134e2056 100644 --- a/GeosCore/fullchem_mod.F90 +++ b/GeosCore/fullchem_mod.F90 @@ -1309,14 +1309,14 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, & ! ! NOTE: KppId is the KPP ID # for each of the prod and loss ! diagnostic species. This is the value used to index the - ! KPP "VAR" array (in module gckpp_Global.F90). + ! KPP "C" array (in module gckpp_Global.F90). !==================================================================== ! Chemical loss of species or families [molec/cm3/s] IF ( State_Diag%Archive_Loss ) THEN DO S = 1, State_Diag%Map_Loss%nSlots KppId = State_Diag%Map_Loss%slot2Id(S) - State_Diag%Loss(I,J,L,S) = VAR(KppID) / DT + State_Diag%Loss(I,J,L,S) = C(KppID) / DT ENDDO ENDIF @@ -1324,7 +1324,7 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, & IF ( State_Diag%Archive_Prod ) THEN DO S = 1, State_Diag%Map_Prod%nSlots KppID = State_Diag%Map_Prod%slot2Id(S) - State_Diag%Prod(I,J,L,S) = VAR(KppID) / DT + State_Diag%Prod(I,J,L,S) = C(KppID) / DT ENDDO ENDIF diff --git a/GeosCore/hco_interface_gc_mod.F90 b/GeosCore/hco_interface_gc_mod.F90 index df4bff748..ee29b21c1 100644 --- a/GeosCore/hco_interface_gc_mod.F90 +++ b/GeosCore/hco_interface_gc_mod.F90 @@ -4627,6 +4627,7 @@ SUBROUTINE Compute_Sflx_for_Vdiff( Input_Opt, State_Chm, State_Diag, & USE HCO_State_GC_Mod, ONLY : ExtState USE HCO_State_GC_Mod, ONLY : HcoState USE Input_Opt_Mod, ONLY : OptInput + USE Mercury_Mod, ONLY : Hg_Emis USE PhysConstants USE Species_Mod, ONLY : Species USE State_Chm_Mod, ONLY : ChmState @@ -4863,6 +4864,11 @@ SUBROUTINE Compute_Sflx_for_Vdiff( Input_Opt, State_Chm, State_Diag, & ENDDO eflx(I,J,NA) = eflx(I,J,NA) + tmpFlx + ! For Hg simulations, also add Hg emissions not handled by HEMCO + IF ( Input_Opt%ITS_A_MERCURY_SIM ) THEN + eflx(I,J,NA) = eflx(I,J,NA) + Hg_EMIS(I,J,NA) + ENDIF + ENDIF #endif @@ -5009,22 +5015,16 @@ SUBROUTINE Compute_Sflx_for_Vdiff( Input_Opt, State_Chm, State_Diag, & IF ( ThisSpc%Is_Hg2 ) THEN - ! Get the category number for this Hg2 tracer - Hg_Cat = ThisSpc%Hg_Cat - ! Archive dry-deposited Hg2 - CALL ADD_Hg2_DD ( I, J, Hg_Cat, dep ) - CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, dep, & + CALL ADD_Hg2_DD ( I, J, dep ) + CALL ADD_Hg2_SNOWPACK( I, J, dep, & State_Met, State_Chm, State_Diag ) ELSE IF ( ThisSpc%Is_HgP ) THEN - ! Get the category number for this HgP tracer - Hg_Cat = ThisSpc%Hg_Cat - ! Archive dry-deposited HgP - CALL ADD_HgP_DD ( I, J, Hg_Cat, dep ) - CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, dep, & + CALL ADD_HgP_DD ( I, J, dep ) + CALL ADD_Hg2_SNOWPACK( I, J, dep, & State_Met, State_Chm, State_Diag ) ENDIF diff --git a/GeosCore/hco_utilities_gc_mod.F90 b/GeosCore/hco_utilities_gc_mod.F90 index 679b7b2d3..73c7377ff 100644 --- a/GeosCore/hco_utilities_gc_mod.F90 +++ b/GeosCore/hco_utilities_gc_mod.F90 @@ -1551,7 +1551,6 @@ SUBROUTINE Get_GC_Restart( Input_Opt, State_Chm, State_Grid, & #ifdef APM USE APM_Init_Mod, ONLY : APMIDS #endif - USE Ocean_Mercury_Mod, ONLY : Check_Ocean_Mercury ! ! !INPUT PARAMETERS: ! @@ -1629,23 +1628,6 @@ SUBROUTINE Get_GC_Restart( Input_Opt, State_Chm, State_Grid, & ! Set pointer to species concentrations Spc => State_Chm%Species - !================================================================= - ! If running Hg simulation, set Hg-specific local variables - !================================================================= - IF ( Input_Opt%ITS_A_MERCURY_SIM ) THEN - - ! Set the # of tagHg categories from State_Chm - Num_Hg_Categories = State_Chm%N_Hg_CATS - - ! Set variable storing names for each of the Hg categories - Hg_Cat_Name => State_Chm%Hg_Cat_Name - - ! Set Hg species index corresponding to a given Hg category number; - ! total is always the first category - Total_Hg_Id = State_Chm%Hg0_Id_List(1) - - ENDIF - !================================================================= ! Open GEOS-Chem restart file !================================================================= @@ -2286,18 +2268,18 @@ SUBROUTINE Get_GC_Restart( Input_Opt, State_Chm, State_Grid, & ! Assign ocean mercury data and write total mass to log file SELECT CASE( M ) - CASE ( 1 ) - State_Chm%OceanHg0(:,:,Total_Hg_Id) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%OceanHg0(:,:,Total_Hg_Id) ), 'kg' - CASE ( 2 ) - State_Chm%OceanHg2(:,:,Total_Hg_Id) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%OceanHg2(:,:,Total_Hg_Id) ), 'kg' - CASE ( 3 ) - State_Chm%OceanHgP(:,:,Total_Hg_Id) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%OceanHgP(:,:,Total_Hg_Id) ), 'kg' + CASE ( 1 ) + State_Chm%OceanHg0 = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%OceanHg0 ), 'kg' + CASE ( 2 ) + State_Chm%OceanHg2 = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%OceanHg2 ), 'kg' + CASE ( 3 ) + State_Chm%OceanHgP = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%OceanHgP ), 'kg' END SELECT ELSE @@ -2309,78 +2291,13 @@ SUBROUTINE Get_GC_Restart( Input_Opt, State_Chm, State_Grid, & ENDDO -!----------------------------------------------------------------------------- -! Tagged Hg simulation has been disabled in 13.4.0 -- the Hg mechanism -! generated via KPP does not use these species any longer: -! -- Bob Yantosca (08 Apr 2022) -! !-------------------------------------------------------------- -! ! Additional tagged ocean Hg species -! !-------------------------------------------------------------- -! IF ( Input_Opt%LSPLIT ) THEN -! DO M = 1, 3 -! DO N = 2, Num_Hg_Categories -! -! ! Define variable name. Include appended region. -! SELECT CASE( M ) -! CASE ( 1 ) -! HgSpc = 'Hg0' -! CASE ( 2 ) -! HgSpc = 'Hg2' -! CASE ( 3 ) -! HgSpc = 'HgP' -! END SELECT -! v_name = 'OCEAN_' // TRIM( HgSpc ) // & -! '_' // TRIM( Hg_Cat_Name(N) ) -! -! ! Get variable from HEMCO and store in local array -! CALL HCO_GC_GetPtr( Input_Opt, State_Grid, TRIM(v_name), & -! Ptr2D, RC, FOUND=FOUND ) -! -! ! Check if variable is in file -! IF ( FOUND ) THEN -! -! ! Assign ocean mercury data and write total mass to log -! SELECT CASE( M ) -! CASE ( 1 ) -! State_Chm%OceanHg0(:,:,N) = Ptr2D -! WRITE( 6, 240 ) TRIM( v_name ), & -! SUM( State_Chm%OceanHg0(:,:,N) ), 'kg' -! CASE ( 2 ) -! State_Chm%OceanHg2(:,:,N) = Ptr2D -! WRITE( 6, 240 ) TRIM( v_name ), & -! SUM( State_Chm%OceanHg2(:,:,N) ), 'kg' -! CASE ( 3 ) -! State_Chm%OceanHgP(:,:,N) = Ptr2D -! WRITE( 6, 240 ) TRIM( v_name ), & -! SUM( State_Chm%OceanHgP(:,:,N) ), 'kg' -! END SELECT -! -! ELSE -! WRITE( 6, 230 ) TRIM( v_name ) -! ENDIF -! -! ! Nullify pointer -! Ptr2D => NULL() -! -! ENDDO -! ENDDO -! -! ! Make sure tagged & total species sum up -! IF ( Input_Opt%USE_CHECKS ) THEN -! CALL CHECK_OCEAN_MERCURY( State_Chm, State_Grid, & -! 'end of READ_GC_RESTART' ) -! ENDIF -! ENDIF -!----------------------------------------------------------------------------- - !-------------------------------------------------------------- ! Hg snowpack on land and ocean !-------------------------------------------------------------- DO M = 1, 4 - DO N = 1, Num_Hg_Categories - ! Define variable name prefix - SELECT CASE( M ) + ! Define variable name prefix + SELECT CASE( M ) CASE ( 1 ) Prefix = 'SNOW_HG_OCEAN' ! Reducible on ocean CASE ( 2 ) @@ -2389,51 +2306,41 @@ SUBROUTINE Get_GC_Restart( Input_Opt, State_Chm, State_Grid, & Prefix = 'SNOW_HG_LAND' ! Reducible on land CASE ( 4 ) Prefix = 'SNOW_HG_LAND_STORED' ! Non-reducible on land - END SELECT + END SELECT - IF ( N == 1 ) THEN - v_name = TRIM( Prefix ) - ELSE - ! Append category name if tagged - v_name = TRIM( Prefix ) // '_' // & - TRIM( Hg_Cat_Name(N) ) - ENDIF + v_name = TRIM( Prefix ) - ! Get variable from HEMCO and store in local array - CALL HCO_GC_GetPtr( Input_Opt, State_Grid, TRIM( v_name ), & - Ptr2D, RC, FOUND=FOUND ) + ! Get variable from HEMCO and store in local array + CALL HCO_GC_GetPtr( Input_Opt, State_Grid, TRIM( v_name ), & + Ptr2D, RC, FOUND=FOUND ) - ! Check if variable is in file - IF ( FOUND ) THEN + ! Check if variable is in file + IF ( FOUND ) THEN - ! Assign ocean mercury data and write total mass to file - SELECT CASE( M ) + ! Assign ocean mercury data and write total mass to file + SELECT CASE( M ) CASE ( 1 ) - State_Chm%SnowHgOcean(:,:,N) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%SnowHgOcean(:,:,N) ), 'kg' + State_Chm%SnowHgOcean = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%SnowHgOcean ), 'kg' CASE ( 2 ) - State_Chm%SnowHgOceanStored(:,:,N) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%SnowHgOceanStored(:,:,N) ),'kg' + State_Chm%SnowHgOceanStored = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%SnowHgOceanStored ),'kg' CASE ( 3 ) - State_Chm%SnowHgLand(:,:,N) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%SnowHgLand(:,:,N) ), 'kg' + State_Chm%SnowHgLand = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%SnowHgLand ), 'kg' CASE ( 4 ) - State_Chm%SnowHgLandStored(:,:,N) = Ptr2D - WRITE( 6, 240 ) TRIM( v_name ), & - SUM( State_Chm%SnowHgLandStored(:,:,N) ), 'kg' + State_Chm%SnowHgLandStored = Ptr2D + WRITE( 6, 240 ) TRIM( v_name ), & + SUM( State_Chm%SnowHgLandStored ), 'kg' END SELECT - ELSE - WRITE( 6, 230 ) TRIM( v_name ) - ENDIF - - ! Nullify pointer - Ptr2D => NULL() + ELSE + WRITE( 6, 230 ) TRIM( v_name ) + ENDIF - ENDDO ENDDO ! Format strings diff --git a/GeosCore/land_mercury_mod.F90 b/GeosCore/land_mercury_mod.F90 index c86075047..57e7aed9b 100644 --- a/GeosCore/land_mercury_mod.F90 +++ b/GeosCore/land_mercury_mod.F90 @@ -44,9 +44,6 @@ MODULE LAND_MERCURY_MOD ! ! !PRIVATE TYPES: ! - ! Keep a shadow variable of N_HG_CATS for backwards compatibility - INTEGER :: N_HG_CATS - ! Plant transpiration rate [m/s] REAL(fp), ALLOCATABLE :: TRANSP(:,:) @@ -83,7 +80,7 @@ SUBROUTINE LAND_MERCURY_FLUX( LFLUX, LHGSNOW, State_Grid, State_Met ) ! !OUTPUT PARAMETERS: ! ! Hg0 flux [kg/s] - REAL(fp), INTENT(OUT) :: LFLUX(State_Grid%NX,State_Grid%NY,N_Hg_CATS) + REAL(fp), INTENT(OUT) :: LFLUX(State_Grid%NX,State_Grid%NY) ! ! !REVISION HISTORY: ! 30 Aug 2010 - N. E. Selin, C. Holmes, B. Corbitt - Initial version @@ -98,7 +95,7 @@ SUBROUTINE LAND_MERCURY_FLUX( LFLUX, LHGSNOW, State_Grid, State_Met ) REAL(fp) :: FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND LOGICAL :: IS_LAND_OR_ICE REAL(fp), PARAMETER :: SEC_PER_YR = 365.25e+0_fp * 86400e+0_fp - INTEGER :: I, J, NN + INTEGER :: I, J !================================================================= ! LAND_MERCURY_FLUX begins here! @@ -107,14 +104,13 @@ SUBROUTINE LAND_MERCURY_FLUX( LFLUX, LHGSNOW, State_Grid, State_Met ) ! Emission timestep [s] DTSRCE = GET_TS_EMIS() - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, NN ) & - !$OMP PRIVATE( REEMFRAC, FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND ) & - !$OMP PRIVATE( IS_LAND_OR_ICE ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, REEMFRAC )& + !$OMP PRIVATE( FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND, IS_LAND_OR_ICE )& + !$OMP COLLAPSE( 2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX - DO NN = 1, N_Hg_CATS ! Distinguish between ice/snow and snow-free land FRAC_SNOW_OR_ICE = MIN( State_Met%FRSNO(I,J) + & @@ -126,6 +122,9 @@ SUBROUTINE LAND_MERCURY_FLUX( LFLUX, LHGSNOW, State_Grid, State_Met ) IS_LAND_OR_ICE = (( FRAC_SNOWFREE_LAND > 0e+0_fp ) .OR. & ( FRAC_SNOW_OR_ICE > 0e+0_fp )) + ! No flux from non-land surfaces (water, sea ice) + LFLUX(I,J) = 0e+0_fp + ! If snow or ice on the ground, reemission fraction is 0.6, ! otherwise 0.2 IF ( IS_LAND_OR_ICE ) THEN @@ -138,20 +137,13 @@ SUBROUTINE LAND_MERCURY_FLUX( LFLUX, LHGSNOW, State_Grid, State_Met ) ENDIF ! Mass of emitted Hg(0), kg - LFLUX(I,J,NN) = ( WD_HgP(I,J,NN) + WD_Hg2(I,J,NN) + & - DD_HgP(I,J,NN) + DD_Hg2(I,J,NN) ) * REEMFRAC + LFLUX(I,J) = ( WD_HgP(I,J) + WD_Hg2(I,J) + & + DD_HgP(I,J) + DD_Hg2(I,J) ) * REEMFRAC ! Emission rate of Hg(0). Convert kg /timestep -> kg/s - LFLUX(I,J,NN) = LFLUX(I,J,NN) / DTSRCE - - ELSE - - ! No flux from non-land surfaces (water, sea ice) - LFLUX(I,J,NN) = 0e+0_fp - + LFLUX(I,J) = LFLUX(I,J) / DTSRCE ENDIF - ENDDO ENDDO ENDDO !$OMP END PARALLEL DO @@ -483,11 +475,12 @@ SUBROUTINE SOILEMIS( EHg0_dist, EHg0_so, State_Grid, State_Met ) ! SOILEMIS begins here! !================================================================= - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, SOIL_EMIS ) & - !$OMP PRIVATE( DRYSOIL_HG, TAUZ, LIGHTFRAC, AREA_M2, SUNCOSVALUE ) & - !$OMP PRIVATE( IS_SNOWFREE_LAND, FRAC_SNOWFREE_LAND ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, SOIL_EMIS )& + !$OMP PRIVATE( DRYSOIL_HG, TAUZ, LIGHTFRAC, AREA_M2, SUNCOSVALUE )& + !$OMP PRIVATE( IS_SNOWFREE_LAND, FRAC_SNOWFREE_LAND )& + !$OMP COLLAPSE( 2 ) DO J=1, State_Grid%NY DO I=1, State_Grid%NX @@ -495,6 +488,9 @@ SUBROUTINE SOILEMIS( EHg0_dist, EHg0_so, State_Grid, State_Met ) State_Met%FRSNO(I,J), 0e+0_fp ) IS_SNOWFREE_LAND = ( FRAC_SNOWFREE_LAND > 0e+0_fp ) + ! First, assume no soil emissions (e.g. for snow/water/ice) + EHg0_so(I,J) = 0e+0_fp + IF ( IS_SNOWFREE_LAND ) THEN ! attenuate solar radiation based on function of leaf area index @@ -526,11 +522,6 @@ SUBROUTINE SOILEMIS( EHg0_dist, EHg0_so, State_Grid, State_Met ) ! Multiply by fractional land area EHg0_so(I,J) = EHg0_so(I,J) * FRAC_SNOWFREE_LAND - ELSE - - ! no soil emissions from water and ice - EHg0_so(I,J) = 0e+0_fp - ENDIF ENDDO @@ -575,7 +566,7 @@ SUBROUTINE SNOWPACK_MERCURY_FLUX( FLUX, LHGSNOW, State_Chm, & ! !OUTPUT PARAMETERS: ! ! Hg0 flux [kg/s] - REAL(fp), INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY,N_Hg_CATS) + REAL(fp), INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY) ! ! !REMARKS: ! Emissions are a linear function of Hg mass stored in the snowpack. The @@ -601,7 +592,7 @@ SUBROUTINE SNOWPACK_MERCURY_FLUX( FLUX, LHGSNOW, State_Chm, & ! ! !LOCAL VARIABLES: ! - INTEGER :: I, J, NN, MONTH + INTEGER :: I, J, MONTH REAL(fp) :: DTSRCE, SNOW_HG_OC_NEW, K_EMIT, SWRAD REAL(fp) :: SNOW_HG_LN_NEW, FLUX_TMP @@ -612,8 +603,8 @@ SUBROUTINE SNOWPACK_MERCURY_FLUX( FLUX, LHGSNOW, State_Chm, & REAL(fp), PARAMETER :: K_EMIT_0 = 2.5e-9_fp ! Pointers - REAL(fp), POINTER :: SNOW_HG_OC(:,:,:) - REAL(fp), POINTER :: SNOW_HG_LN(:,:,:) + REAL(fp), POINTER :: SNOW_HG_OC(:,:) + REAL(fp), POINTER :: SNOW_HG_LN(:,:) !================================================================= ! SNOWPACK_MERCURY_FLUX begins here! @@ -632,11 +623,11 @@ SUBROUTINE SNOWPACK_MERCURY_FLUX( FLUX, LHGSNOW, State_Chm, & ! Emission timestep [s] DTSRCE = GET_TS_EMIS() - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, NN ) & - !$OMP PRIVATE( SNOW_HG_OC_NEW, K_EMIT, SWRAD, SNOW_HG_LN_NEW ) & - !$OMP PRIVATE( FLUX_TMP ) + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, SNOW_HG_OC_NEW, K_EMIT )& + !$OMP PRIVATE( SWRAD, SNOW_HG_LN_NEW, FLUX_TMP )& + !$OMP COLLAPSE( 2 ) DO J = 1, State_Grid%NY DO I = 1, State_Grid%NX @@ -659,40 +650,36 @@ SUBROUTINE SNOWPACK_MERCURY_FLUX( FLUX, LHGSNOW, State_Chm, & ! K_EMIT if it is >= K_EMIT_0 (jaf, 5/19/11) IF ( K_EMIT >= K_EMIT_0 ) THEN - ! Loop over # of categories - DO NN = 1, N_Hg_CATS + ! Zero temporary reservoir on each iteration + FLUX_TMP = 0D0 - ! Zero temporary reservoir on each iteration - FLUX_TMP = 0D0 + ! Check if there is Hg that could be emitted + IF ( SNOW_HG_OC(I,J) > 0e+0_fp ) THEN - ! Check if there is Hg that could be emitted - IF ( SNOW_HG_OC(I,J,NN) > 0e+0_fp ) THEN + ! New mass of snow in Hg + SNOW_HG_OC_NEW = SNOW_HG_OC(I,J) * EXP(-K_EMIT*DTSRCE) - ! New mass of snow in Hg - SNOW_HG_OC_NEW = SNOW_HG_OC(I,J,NN) * EXP(-K_EMIT*DTSRCE) + FLUX_TMP = MAX(SNOW_HG_OC(I,J) - SNOW_HG_OC_NEW, 0e+0_fp) - FLUX_TMP = MAX(SNOW_HG_OC(I,J,NN) - SNOW_HG_OC_NEW, 0e+0_fp) + SNOW_HG_OC(I,J) = SNOW_HG_OC_NEW - SNOW_HG_OC(I,J,NN) = SNOW_HG_OC_NEW + ENDIF !SNOW_HG_OC > 0 + + ! Check if there is Hg in land snow that could be emitted + IF ( SNOW_HG_LN(I,J) > 0e+0_fp ) THEN - ENDIF !SNOW_HG_OC > 0 + ! New mass of snow in Hg + SNOW_HG_LN_NEW = SNOW_HG_LN(I,J) * EXP(-K_EMIT*DTSRCE) - ! Check if there is Hg in land snow that could be emitted - IF ( SNOW_HG_LN(I,J,NN) > 0e+0_fp ) THEN + FLUX_TMP = FLUX_TMP + MAX(SNOW_HG_LN(I,J)-SNOW_HG_LN_NEW,0D0) - ! New mass of snow in Hg - SNOW_HG_LN_NEW = SNOW_HG_LN(I,J,NN) * EXP(-K_EMIT*DTSRCE) + SNOW_HG_LN(I,J) = SNOW_HG_LN_NEW - FLUX_TMP = FLUX_TMP + MAX(SNOW_HG_LN(I,J,NN)-SNOW_HG_LN_NEW,0D0) + ENDIF !SNOW_HG_LN > 0 - SNOW_HG_LN(I,J,NN) = SNOW_HG_LN_NEW + ! Convert mass -> flux + FLUX(I,J) = FLUX_TMP / DTSRCE - ENDIF !SNOW_HG_LN > 0 - - ! Convert mass -> flux - FLUX(I,J,NN) = FLUX_TMP / DTSRCE - - ENDDO !Loop over NN ENDIF !K_EMIT >= K_EMIT_0 ENDDO @@ -1004,9 +991,6 @@ SUBROUTINE INIT_LAND_MERCURY( Input_Opt, State_Chm, RC ) ! Assume success RC = GC_SUCCESS - ! Save State_Chm%N_Hg_CATS in a shadow variable - N_Hg_CATS = State_Chm%N_Hg_CATS - END SUBROUTINE INIT_LAND_MERCURY !EOC !------------------------------------------------------------------------------ diff --git a/GeosCore/mercury_mod.F90 b/GeosCore/mercury_mod.F90 index 0ff39b149..ecf9b43db 100644 --- a/GeosCore/mercury_mod.F90 +++ b/GeosCore/mercury_mod.F90 @@ -48,6 +48,7 @@ MODULE Mercury_Mod ! PUBLIC :: Cleanup_Mercury PUBLIC :: ChemMercury + PUBLIC :: EmissMercury PUBLIC :: Init_Mercury ! ! !REMARKS: @@ -130,6 +131,11 @@ MODULE Mercury_Mod INTEGER :: Map_Hg2gas(25) INTEGER, ALLOCATABLE :: PL_Kpp_ID(:) REAL(fp), ALLOCATABLE :: EHg0_an(:,:) + REAL(fp), ALLOCATABLE :: EHg0_dist(:,:) + REAL(fp), ALLOCATABLE :: EHg0_ln(:,:) + REAL(fp), ALLOCATABLE :: EHg0_oc(:,:) + REAL(fp), ALLOCATABLE :: EHg0_so(:,:) + REAL(fp), ALLOCATABLE :: EHg0_snow(:,:) REAL(fp), ALLOCATABLE :: EHg2_an(:,:) REAL(fp), ALLOCATABLE :: COSZM(:,:) ! Max daily SZA REAL(fp), ALLOCATABLE :: srMw(:) @@ -140,6 +146,11 @@ MODULE Mercury_Mod CHARACTER(LEN=8), & ALLOCATABLE :: AerSpcNames(:) + ! For now, we need an emission array for the HG simulation + ! that can be passed to vdiff_mod.F90 since Trac_Tend does + ! not exist anymore (ckeller, 10/21/2014). + REAL(fp), ALLOCATABLE, PUBLIC :: HG_EMIS(:,:,:) + !-------------------------------------------------------------------------- ! Pointers to fields in the HEMCO data structure, which must be REAL*4. ! (NOTE: We can set them to NULL here because hey are globally @@ -189,6 +200,401 @@ MODULE Mercury_Mod !------------------------------------------------------------------------------ !BOP ! +! !IROUTINE: emissmercury +! +! !DESCRIPTION: Subroutine EMISSMERCURY is the driver routine for mercury +! emissions. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE EmissMercury( Input_Opt, State_Chm, State_Diag, & + State_Grid, State_Met, RC ) +! +! !USES: +! + USE DEPO_MERCURY_MOD, ONLY : RESET_HG_DEP_ARRAYS + USE ErrCode_Mod + USE ERROR_MOD + USE HCO_Utilities_GC_Mod, ONLY : HCO_GC_GetPtr + USE Input_Opt_Mod, ONLY : OptInput + USE Land_Mercury_Mod, ONLY : Land_Mercury_Flux + USE Land_Mercury_Mod, ONLY : SOILEMIS + USE Land_Mercury_Mod, ONLY : SNOWPACK_MERCURY_FLUX + USE Ocean_Mercury_Mod, ONLY : OCEAN_MERCURY_FLUX + USE State_Chm_Mod, ONLY : ChmState + USE State_Diag_Mod, ONLY : DgnState + USE State_Grid_Mod, ONLY : GrdState + USE State_Met_Mod, ONLY : MetState + USE Time_Mod, ONLY : GET_MONTH, ITS_A_NEW_MONTH + USE UnitConv_Mod, ONLY : Convert_Spc_Units + + ! Added for GTMM (ccc, 11/19/09) + !USE LAND_MERCURY_MOD, ONLY : GTMM_DR +! +! !INPUT PARAMETERS: +! + TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object + TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object + TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object +! +! !INPUT/OUTPUT PARAMETERS: +! + TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object + TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diagnostics State object +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT(OUT) :: RC ! Success or failure? +! +! !REMARKS: +! +! +! !REVISION HISTORY: +! 03 Jun 2013 - N. (Eckley) Selin - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + LOGICAL, SAVE :: FIRST = .TRUE. + INTEGER :: THISMONTH, I, J + LOGICAL :: prtDebug + + ! Pointers + REAL(f4), POINTER :: Ptr2D(:,:) + + ! Strings + CHARACTER(LEN=63) :: OrigUnit + CHARACTER(LEN=255) :: ErrMsg + CHARACTER(LEN=255) :: ThisLoc + + !================================================================= + ! EMISSMERCURY begins here! + !================================================================= + + ! Assume success + RC = GC_SUCCESS + prtDebug = ( Input_Opt%LPRT .and. Input_Opt%amIRoot ) + ErrMsg = '' + ThisLoc = ' -> at EMISSMERCURY (in module GeosCore/mercury_mod.F90)' + + ! Convert species units to [kg] for EMISSMERCURY (ewl, 8/12/15) + CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, & + 'kg', RC, OrigUnit=OrigUnit ) + + ! Trap potential errors + IF ( RC /= GC_SUCCESS ) THEN + ErrMsg = 'Error encountered in "Convert_Spc_Units" #1!' + CALL GC_Error( ErrMsg, RC, ThisLoc ) + RETURN + ENDIF + + !================================================================= + ! Get data pointers from HEMCO on the first call + !================================================================= + IF ( FIRST ) THEN + + ! Soil distribution + CALL HCO_GC_GetPtr( Input_Opt, State_Grid, 'HG0_SOILDIST', Ptr2D, RC ) + IF ( RC /= GC_SUCCESS ) THEN + ErrMsg = 'Could not get pointer to HEMCO field HG0_SOILDIST!' + CALL GC_Error( ErrMsg, RC, ThisLoc ) + RETURN + ENDIF + EHg0_dist = Ptr2D(:,:) + Ptr2D => NULL() + + ! Reset first-time flag + FIRST = .FALSE. + ENDIF + + !-------------------------------------------------------------- + ! Here we are NOT using the Global Terrstrial Mercury Model... + !-------------------------------------------------------------- + + ! Read offline ocean Hg0 + EHg0_oc = 0.0_fp + CALL OFFLINEOCEAN_READMO( State_Chm, State_Diag, State_Grid, & + State_Met, EHg0_oc, RC ) + + ! Get land mercury flux of Hg0 + CALL LAND_MERCURY_FLUX( LHgSNOW = LHgSNOW, & + State_Grid = State_Grid, & + State_Met = State_Met, & + LFLUX = EHg0_ln ) + IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a LAND_FLUX' ) + + ! Get soil mercury flux of Hg0 + CALL SOILEMIS( EHg0_dist = EHg0_dist, & + State_Grid = State_Grid, & + State_Met = State_Met, & + EHg0_so = EHg0_so ) + IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SOILEMIS' ) + + ! Get snow mercury flux of Hg0 + CALL SNOWPACK_MERCURY_FLUX( LHgSNOW = LHgSNOW, & + State_Chm = State_Chm, & + State_Grid = State_Grid, & + State_Met = State_Met, & + FLUX = EHg0_snow ) + IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SNOW_FLUX' ) + + ! If we are using the non-local PBL mixing, + ! we need to initialize the EMIS_SAVE array (cdh, 08/27/09) + ! EMIS_SAVE is now HG_EMIS (ckeller, 10/21/2014) + IF ( Input_Opt%LNLPBL ) HG_EMIS = 0.0e+0_fp + + ! Zero arrays for Hg deposition + CALL RESET_HG_DEP_ARRAYS() + IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: ' // & + 'a RESET_HG_DEP_ARRAYS' ) + + ! Add Hg(0) source into State_Chm%Species [kg] + CALL SRCHg0( Input_Opt, State_Chm, State_Diag, State_Grid, State_Met, RC ) + IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SRCHg0' ) + + ! Convert species units back to original unit + CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, & + OrigUnit, RC ) + IF ( RC /= GC_SUCCESS ) THEN + CALL GC_Error('Unit conversion error', RC, & + 'Routine EMISSMERCURY in mercury_mod.F90') + RETURN + ENDIF + + END SUBROUTINE EMISSMERCURY +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: emithg +! +! !DESCRIPTION: Subroutine EMITHG directs emission either to the chemical +! species array (State\_Chm%Species) directly or to Hg\_EMIS for use by the +! non-local PBL mixing. This is a programming convenience. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE EMITHG( I, J, L, N, E_HG, Input_Opt, State_Chm, State_Grid ) +! +! !USES: +! + USE ErrCode_Mod + USE Input_Opt_Mod, ONLY : OptInput + USE Species_Mod, ONLY : SpcConc + USE State_Chm_Mod, ONLY : ChmState + USE State_Grid_Mod, ONLY : GrdState + USE TIME_MOD, ONLY : GET_TS_EMIS +! +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN) :: I, J, L, N ! Grid boxes + species # + REAL(fp), INTENT(IN) :: E_Hg ! Hg emissions + TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object + TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object +! +! !INPUT/OUTPUT PARAMETERS: +! + TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object +! +! !REVISION HISTORY: +! 27 Aug 2009 - C. Holmes - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + ! Scalars + REAL(fp) :: AM2, TS + + ! Pointers + TYPE(SpcConc), POINTER :: Spc(:) + + !================================================================= + ! EMITHG begins here! + !================================================================= + + IF ( Input_Opt%LNLPBL ) THEN + + !-------------------------------------------------------------- + ! We are using FULL PBL MIXING (routine TURBDAY) + ! + ! Save emissions for non-local PBL mixing or emit directly. + ! Make sure that emitted mass is non-negative + ! This is hear only for consistency with old code which warned + ! of underflow error (cdh, 08/27/09) + ! EMIS_SAVE is now HG_EMIS array. Convert kg to kg/m2/s + ! (ckeller, 09/23/2014) + !-------------------------------------------------------------- + + ! Surface area [m2] + AM2 = State_Grid%Area_M2(I,J) + + ! Emission timestep + TS = GET_TS_EMIS() + + ! Save emissions as [kg/m2/s]. These will be added + ! to the chemical species array in routine DO_TEND + ! (in mixing_mod.F90). + Hg_EMIS(I,J,N) = Hg_EMIS(I,J,N) + ( MAX( E_Hg, 0.0_fp ) / AM2 / TS ) + + ELSE + + !-------------------------------------------------------------- + ! We are using FULL PBL MIXING (routine TURBDAY) + ! so add directly to the State_Chm%Species array + !-------------------------------------------------------------- + + ! Point to the chemical spcies array [kg] + Spc => State_Chm%Species + + ! Add emissions + Spc(N)%Conc(I,J,L) = Spc(N)%Conc(I,J,L) + MAX( E_Hg, 0.0_fp ) + + ! Free pointer + Spc => NULL() + + ENDIF + + END SUBROUTINE EMITHG +!EOP +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: srcHg0 +! +! !DESCRIPTION: Subroutine SRCHg0 is the subroutine for Hg(0) emissions. +! Emissions of Hg(0) will be distributed throughout the boundary layer. +! (eck, cdh, bmy, 1/21/05, 4/6/06) +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE SRCHg0( Input_Opt, State_Chm, State_Diag, & + State_Grid, State_Met, RC ) +! +! !USES: +! + USE ErrCode_Mod + USE Input_Opt_Mod, ONLY : OptInput + USE State_Chm_Mod, ONLY : ChmState + USE State_Diag_Mod, ONLY : DgnState + USE State_Grid_Mod, ONLY : GrdState + USE State_Met_Mod, ONLY : MetState + USE TIME_MOD, ONLY : GET_TS_EMIS +! +! !INPUT PARAMETERS: +! + TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object + TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object + TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object +! +! !INPUT/OUTPUT PARAMETERS: +! + TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object + TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diagnostics State object +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT(OUT) :: RC ! Success or failure? +! +! !REVISION HISTORY: +! 21 Jan 2005 - N. (Eckley) Selin, C. Holmes - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + INTEGER :: I, J, L + REAL(fp) :: DTSRCE, E_Hg, T_Hg + + !================================================================= + ! SRCHg0 begins here! + !================================================================= + + ! Initialize + RC = GC_SUCCESS + DTSRCE = GET_TS_EMIS() ! Timestep [s] + + ! Zero ocean and snow if the "anthro Hg only" option is selected + IF ( LAnthroHgOnly ) THEN + EHg0_oc = 0.0_fp + EHg0_snow = 0.0_fp + ENDIF + + ! Loop over grid boxes + !$OMP PARALLEL DO & + !$OMP DEFAULT( SHARED )& + !$OMP PRIVATE( I, J, L, T_Hg, E_Hg )& + !$OMP COLLAPSE( 2 ) + DO J = 1, State_Grid%NY + DO I = 1, State_Grid%NX + + ! Compute total Hg(0) emissions not computed by HEMCO + ! (oceans + land + natural + snow) + T_Hg = EHg0_oc(I,J) + EHg0_ln(I,J) + EHg0_so(I,J) + EHg0_snow(I,J) + + !===================================================================== + ! Partition Hg0 throughout PBL; store into State_Chm%Species [kg] + ! Now make sure State_Chm%Species does not underflow (cdh, bmy, 4/6/06) + !===================================================================== + DO L = 1, State_Met%PBL_MAX_L + E_Hg = State_Met%F_OF_PBL(I,J,L) * T_Hg * DTSRCE + CALL EmitHg( I, J, L, id_Hg0, E_Hg, Input_Opt, State_Chm, State_Grid ) + ENDDO + + !===================================================================== + ! %%%%% HISTORY (aka netCDF diagnostics) %%%%% + ! + ! Save the various Hg0 emissions (land, ocean, snow, soil) + ! in units of [kg/s]. + ! + ! NOTE: All other emission categories are archived via HEMCO. + !===================================================================== + + ! Land Hg0 emissions [kg/s] + IF ( State_Diag%Archive_EmisHg0land ) THEN + State_Diag%EmisHg0land(I,J) = EHg0_ln(I,J) + ENDIF + + ! Oceanic Hg0 emissions [kg/s] + IF ( State_Diag%Archive_EmisHg0ocean ) THEN + State_Diag%EmisHg0ocean(I,J) = EHg0_oc(I,J) + ENDIF + + ! Snow Hg0 emissions [kg/s] + IF ( State_Diag%Archive_EmisHg0snow ) THEN + State_Diag%EmisHg0snow(I,J) = EHg0_snow(I,J) + ENDIF + + ! Soil Hg0 emissions [kg/s] + IF ( State_Diag%Archive_EmisHg0soil ) THEN + State_Diag%EmisHg0soil(I,J) = EHg0_so(I,J) + ENDIF + + ENDDO + ENDDO + !$OMP END PARALLEL DO + + END SUBROUTINE SRCHg0 +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! ! !IROUTINE: chemmercury ! ! !DESCRIPTION: Subroutine CHEMMERCURY is the driver routine for mercury @@ -2511,8 +2917,8 @@ SUBROUTINE PartitionHg2( Input_Opt, State_Chm, State_Diag, & !================================================================== ! Concentration of Hg2 on aerosols [molec/cm3] - aerConc = Spc(id_Hg2STRP)%Conc(I,J,L) & - + Spc(id_Hg2ORGP)%Conc(I,J,L) & + aerConc = Spc(id_Hg2STRP)%Conc(I,J,L) & + + Spc(id_Hg2ORGP)%Conc(I,J,L) & + Spc(id_Hg2ClP)%Conc(I,J,L) ! Zero organic Hg2 aerosol and Hg2Cl aerosol in stratosphere @@ -2756,7 +3162,7 @@ SUBROUTINE OfflineOcean_ReadMo( State_Chm, State_Diag, State_Grid, & ! ! !OUTPUT PARAMETERS: ! - REAL*8, INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY,1) + REAL*8, INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY) INTEGER, INTENT(OUT) :: RC ! Success or failure? ! ! !REVISION HISTORY: @@ -2780,9 +3186,9 @@ SUBROUTINE OfflineOcean_ReadMo( State_Chm, State_Diag, State_Grid, & REAL(fp) :: TC, TK, Kw REAL(fp) :: Sc, ScCO2, USQ REAL(fp) :: FRAC_L, FRAC_O, H, D - REAL(fp) :: FUP(State_Grid%NX,State_Grid%NY,N_Hg_CATS) - REAL(fp) :: FDOWN(State_Grid%NX,State_Grid%NY,N_Hg_CATS) - REAL(fp) :: Hg0aq(State_Grid%NX,State_Grid%NY,N_Hg_CATS) + REAL(fp) :: FUP(State_Grid%NX,State_Grid%NY) + REAL(fp) :: FDOWN(State_Grid%NX,State_Grid%NY) + REAL(fp) :: Hg0aq(State_Grid%NX,State_Grid%NY) REAL(fp) :: MHg0_air ! Conversion factor from [cm/h * ng/L] --> [kg/m2/s] @@ -2843,22 +3249,17 @@ SUBROUTINE OfflineOcean_ReadMo( State_Chm, State_Diag, State_Grid, & ENDIF ! Only doing Hg0 overall, should add trap for LSPLIT (cpt - 2017) - DO NNN = 1, N_Hg_CATS - IF ( NNN .eq. 1 .and. associated(OCEAN_CONC)) THEN - Hg0aq(:,:,NNN) = OCEAN_CONC(:,:,1) - ELSE - Hg0aq(:,:,NNN) = 0.0e+0_fp - ENDIF - ENDDO + Hg0aq = OCEAN_CONC(:,:,1) ! Emission timestep [s] DTSRCE = GET_TS_EMIS() ! Loop over surface boxes + ! NOTE: Remove the loop over NN -- the tagged Hg simulation is not used !$OMP PARALLEL DO & !$OMP DEFAULT( SHARED )& !$OMP PRIVATE( I, A_M2, vi, ScCO2 )& - !$OMP PRIVATE( J, NN, TC, TK )& + !$OMP PRIVATE( J, TC, TK )& !$OMP PRIVATE( N, CHg0, FRAC_L, FRAC_O )& !$OMP PRIVATE( H, Kw, CHg0aq, Hg0aqtemp, MHg0_air )& !$OMP PRIVATE( IS_OCEAN_BOX, Sc, Usq, D )& @@ -2938,96 +3339,85 @@ SUBROUTINE OfflineOcean_ReadMo( State_Chm, State_Diag, State_Grid, & ! Mass transfer coefficient [cm/h], from Nightingale et al. 2000 Kw = ( 0.25_fp * Usq ) / SQRT( Sc / ScCO2 ) - ! Loop over all Hg categories - DO NN = 1, N_Hg_CATS + ! Hg0 tracer number (for Spc) + N = id_Hg0 !Hg0_Id_List(NN) - ! Hg0 tracer number (for Spc) - N = id_Hg0 !Hg0_Id_List(NN) + !-------------------------------------------------------- + ! Calculate oceanic and gas-phase concentration of Hg(0) + !-------------------------------------------------------- - !-------------------------------------------------------- - ! Calculate oceanic and gas-phase concentration of Hg(0) - !-------------------------------------------------------- + ! Concentration of Hg(0) in the ocean [ng/L] + ! now converting from Hg0aq in mol/m3 to ng/L + CHg0aq = Hg0aq(I,J) * 200.59_fp * 1.0e9_fp / 1.0e3_fp - ! Concentration of Hg(0) in the ocean [ng/L] - ! now converting from Hg0aq in mol/m3 to ng/L - CHg0aq = Hg0aq(I,J,NN) * 200.59_fp * 1.0e9_fp / 1.0e3_fp + ! Gas phase Hg(0) concentration: convert [kg] -> [ng/L] + MHg0_air = State_Chm%Species(N)%Conc(I,J,1) + CHg0 = MHg0_air * 1.0e9_fp /State_Met%AIRVOL(I,J,1) - ! Gas phase Hg(0) concentration: convert [kg] -> [ng/L] - MHg0_air = State_Chm%Species(N)%Conc(I,J,1) - CHg0 = MHg0_air * 1.0e9_fp /State_Met%AIRVOL(I,J,1) + !-------------------------------------------------------- + ! Compute flux of Hg(0) from the ocean to the air + !-------------------------------------------------------- - !-------------------------------------------------------- - ! Compute flux of Hg(0) from the ocean to the air - !-------------------------------------------------------- + ! Compute ocean flux of Hg0 [cm/h*ng/L] + FLUX(I,J) = Kw * ( CHg0aq - ( CHg0 / H ) ) - ! Compute ocean flux of Hg0 [cm/h*ng/L] - FLUX(I,J,NN) = Kw * ( CHg0aq - ( CHg0 / H ) ) + !Extra diagnostic: compute flux up and flux down + FUP(I,J) = ( Kw * CHg0aq ) + FDOWN(I,J) = ( Kw * CHg0 / H ) - !Extra diagnostic: compute flux up and flux down - FUP(I,J,NN) = ( Kw * CHg0aq ) - FDOWN(I,J,NN) = ( Kw * CHg0 / H ) + !-------------------------------------------------- + ! Convert [cm/h*ng/L] --> [kg/m2/s] --> [kg/s] + ! Also account for ocean fraction of grid box + FLUX(I,J) = FLUX(I,J) * TO_KGM2S * A_M2 * FRAC_O - !-------------------------------------------------- - ! Convert [cm/h*ng/L] --> [kg/m2/s] --> [kg/s] - ! Also account for ocean fraction of grid box - FLUX(I,J,NN) = FLUX(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O + ! hmh 5/11/16 reverting to old version and uncommenting here + FUP(I,J) = FUP(I,J) * TO_KGM2S * A_M2 * FRAC_O + FDOWN(I,J) = FDOWN(I,J) * TO_KGM2S * A_M2 * FRAC_O + !-------------------------------------------------- - ! hmh 5/11/16 reverting to old version and uncommenting here - FUP(I,J,NN) = FUP(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O - FDOWN(I,J,NN) = FDOWN(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O - !-------------------------------------------------- + !-------------------------------------------------------- + ! Flux limited by ocean and atm Hg(0) + !-------------------------------------------------------- - !-------------------------------------------------------- - ! Flux limited by ocean and atm Hg(0) - !-------------------------------------------------------- + ! Cap the flux w/ the available Hg(0) ocean mass + Hg0aqtemp = CHg0aq * A_M2 * FRAC_O *1.0e-8_fp - ! Cap the flux w/ the available Hg(0) ocean mass - ! THIS IS THE PROBLEM! - Hg0aqtemp = CHg0aq * A_M2 * FRAC_O *1.0e-8_fp - - IF ( FLUX(I,J,NN) * DTSRCE > Hg0aqtemp ) THEN - FLUX(I,J,NN) = Hg0aqtemp / DTSRCE - FUP(I,J,NN) = FLUX(I,J,NN)-FDOWN(I,J,NN) - ENDIF - - ! Cap the neg flux w/ the available Hg(0) atm mass - IF ( (-FLUX(I,J,NN) * DTSRCE ) > MHg0_air ) THEN - FLUX(I,J,NN) = -MHg0_air / DTSRCE - ENDIF - - ! make sure Fup and Fdown do not underflow either - ! debug 2x2.5 diagnostic? - FUP(I,J,NN) = MAX(FUP(I,J,NN), SMALLNUM ) - FDOWN(I,J,NN) = MAX(FDOWN(I,J,NN), SMALLNUM ) + IF ( FLUX(I,J) * DTSRCE > Hg0aqtemp ) THEN + FLUX(I,J) = Hg0aqtemp / DTSRCE + FUP(I,J) = FLUX(I,J)-FDOWN(I,J) + ENDIF - !-------------------------------------------------------- - ! %%%%% HISTORY (aka netCDF diagnostics) %%%%% - ! - ! Fluxes of Hg0 from air to ocean and ocean to air [kg/s] - ! NOTE: Implement for total Hg species at ths time - !-------------------------------------------------------- - IF ( NN == 1 ) THEN + ! Cap the neg flux w/ the available Hg(0) atm mass + IF ( (-FLUX(I,J) * DTSRCE ) > MHg0_air ) THEN + FLUX(I,J) = -MHg0_air / DTSRCE + ENDIF - ! Flux of Hg0 from ocean to air [kg/s] - IF ( State_Diag%Archive_FluxHg0fromOceanToAir ) THEN - State_Diag%FluxHg0fromOceanToAir(I,J) = FUP(I,J,NN) - ENDIF + ! make sure Fup and Fdown do not underflow either + ! debug 2x2.5 diagnostic? + FUP(I,J) = MAX( FUP(I,J), SMALLNUM ) + FDOWN(I,J) = MAX( FDOWN(I,J), SMALLNUM ) - IF ( State_Diag%Archive_FluxHg0fromAirToOcean ) THEN - State_Diag%FluxHg0fromAirToOcean(I,J) = FDOWN(I,J,NN) - ENDIF + !-------------------------------------------------------- + ! %%%%% HISTORY (aka netCDF diagnostics) %%%%% + ! + ! Fluxes of Hg0 from air to ocean and ocean to air [kg/s] + ! NOTE: Implement for total Hg species at ths time + !-------------------------------------------------------- - ENDIF + ! Flux of Hg0 from ocean to air [kg/s] + IF ( State_Diag%Archive_FluxHg0fromOceanToAir ) THEN + State_Diag%FluxHg0fromOceanToAir(I,J) = FUP(I,J) + ENDIF - ENDDO + IF ( State_Diag%Archive_FluxHg0fromAirToOcean ) THEN + State_Diag%FluxHg0fromAirToOcean(I,J) = FDOWN(I,J) + ENDIF ELSE - DO NN = 1, N_Hg_CATS - FLUX(I,J,NN) = 0.0_fp - FUP(I,J,NN) = 0.0_fp - FDOWN(I,J,NN) = 0.0_fp - ENDDO + FLUX(I,J) = 0.0_fp + FUP(I,J) = 0.0_fp + FDOWN(I,J) = 0.0_fp ENDIF ENDDO @@ -3118,7 +3508,6 @@ SUBROUTINE Init_Mercury( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) SpcInfo => NULL() errMsg = '' thisLoc = '-> at DEFINE_TAGGED_Hg (in GeosCore/mercury_mod.F90)' - N_Hg_CATS = 1 ! Write reactions WRITE( 6 ,'(a)' ) ' KPP Reaction Reference ' @@ -3280,6 +3669,31 @@ SUBROUTINE Init_Mercury( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) IF ( RC /= GC_SUCCESS ) RETURN EHg0_an = 0e+0_fp + ALLOCATE( EHg0_dist( State_Grid%NX, State_Grid%NY), STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:EHg0_dist', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + EHg0_dist = 0e+0_fp + + ALLOCATE( EHg0_ln( State_Grid%NX, State_Grid%NY ), STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:EHg0_ln', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + EHg0_ln = 0e+0_fp + + ALLOCATE( EHg0_oc( State_Grid%NX, State_Grid%NY ), STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:EHg0_oc', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + EHg0_oc = 0e+0_fp + + ALLOCATE( EHg0_so( State_Grid%NX, State_Grid%NY ), STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:EHg0_snow', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + EHg0_so = 0e+0_fp + + ALLOCATE( EHg0_snow( State_Grid%NX, State_Grid%NY ), STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:EHg0_snow', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + EHg0_snow = 0e+0_fp + ALLOCATE( EHg2_an( State_Grid%NX, State_Grid%NY ), STAT=RC ) CALL GC_CheckVar( 'mercury_mod.F90:EHg2_an', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN @@ -3314,6 +3728,13 @@ SUBROUTINE Init_Mercury( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) IF ( RC /= GC_SUCCESS ) RETURN srMw = 0.0_fp + ! HG_EMIS is needed for non-local PBL mixing + ALLOCATE( HG_EMIS( State_Grid%NX, State_Grid%NY, State_Chm%nAdvect ), & + STAT=RC ) + CALL GC_CheckVar( 'mercury_mod.F90:HG_EMIS', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + HG_EMIS = 0e+0_fp + !======================================================================== ! Allocate and initialize oxidant concentration pointer !======================================================================== @@ -3726,12 +4147,17 @@ SUBROUTINE Cleanup_Mercury ! Deallocate module arrays IF ( ALLOCATED( COSZM ) ) DEALLOCATE( COSZM ) IF ( ALLOCATED( EHg0_an ) ) DEALLOCATE( EHg0_an ) + IF ( ALLOCATED( EHg0_ln ) ) DEALLOCATE( EHg0_ln ) + IF ( ALLOCATED( EHg0_oc ) ) DEALLOCATE( EHg0_oc ) + IF ( ALLOCATED( EHg0_so ) ) DEALLOCATE( EHg0_so ) + IF ( ALLOCATED( EHg0_snow ) ) DEALLOCATE( EHg0_snow ) IF ( ALLOCATED( EHg2_an ) ) DEALLOCATE( EHg2_an ) IF ( ALLOCATED( srMw ) ) DEALLOCATE( srMw ) IF ( ALLOCATED( TCOSZ ) ) DEALLOCATE( TCOSZ ) IF ( ALLOCATED( TTDAY ) ) DEALLOCATE( TTDAY ) IF ( ALLOCATED( ZERO_DVEL ) ) DEALLOCATE( ZERO_DVEL ) IF ( ALLOCATED( HG2_SEASALT_LOSSRATE ) ) DEALLOCATE( HG2_SEASALT_LOSSRATE ) + IF ( ALLOCATED( HG_EMIS ) ) DEALLOCATE( HG_EMIS ) ! Free pointers to HEMCO fields O3 => NULL() diff --git a/GeosCore/ocean_mercury_mod.F90 b/GeosCore/ocean_mercury_mod.F90 index cfc0ff391..6c3f30369 100644 --- a/GeosCore/ocean_mercury_mod.F90 +++ b/GeosCore/ocean_mercury_mod.F90 @@ -27,7 +27,6 @@ MODULE OCEAN_MERCURY_MOD ! !PUBLIC MEMBER FUNCTIONS: ! PUBLIC :: INIT_OCEAN_MERCURY - PUBLIC :: CHECK_OCEAN_MERCURY PUBLIC :: CLEANUP_OCEAN_MERCURY PUBLIC :: OCEAN_MERCURY_FLUX PUBLIC :: LDYNSEASALT, LGCAPEMIS, LPOLARBR, LBRCHEM, LBROCHEM @@ -94,10 +93,6 @@ MODULE OCEAN_MERCURY_MOD ! !============================================================================ ! Hg_RST_FILE : Name of restart file with ocean tracers - ! USE_CHECKS : Flag for turning on error-checking - ! MAX_RELERR : Max error for total-tag error check [unitless] - ! MAX_ABSERR : Max abs error for total-tag err chk [unitless] - ! MAX_FLXERR : Max error tol for flux error check [unitless] ! Hg2aq_tot : Total Hg2 conc. in the mixed layer [kg ] ! DD_Hg2 : Array for Hg(II) dry dep'd to ocean [kg ] ! Hgaq_tot : Total Hg conc. in the mixed layer [kg ] @@ -145,9 +140,9 @@ MODULE OCEAN_MERCURY_MOD REAL(fpp) :: FLOWNOW ! Private arrays - REAL(fpp), ALLOCATABLE :: Hgaq_tot(:,:,:) + REAL(fpp), ALLOCATABLE :: Hgaq_tot(:,:) REAL(fpp), ALLOCATABLE :: dMLD(:,:) - REAL(fpp), ALLOCATABLE :: HgPaq_SUNK(:,:,:) + REAL(fpp), ALLOCATABLE :: HgPaq_SUNK(:,:) REAL(fpp), ALLOCATABLE :: MLDav(:,:) REAL(fpp), ALLOCATABLE :: newMLD(:,:) REAL(fpp), ALLOCATABLE :: prevMLD(:,:) @@ -207,25 +202,6 @@ MODULE OCEAN_MERCURY_MOD REAL(f4), POINTER :: dMLD1 (:,:) => NULL() REAL(f4), POINTER :: dMLD2 (:,:) => NULL() - !----------------------------------------------------------------- - ! Now keep local Hg tracer indexing variables instead of getting - ! them from tracerid_mod.f. The tracerid_mod.f module is going - ! to be removed for the FlexChem implementation. (bmy, 4/26/16) - !----------------------------------------------------------------- - - ! Scalars for Hg indexing - INTEGER :: N_Hg_CATS - INTEGER :: ID_Hg_atl, ID_Hg_nat, ID_Hg_sat - INTEGER :: ID_Hg_npa, ID_Hg_arc, ID_Hg_ant - INTEGER :: ID_Hg_ocn, ID_Hg_tot - - ! Pointers for Hg indexing - ! NOTE: We can declare these NULL here because - ! these are SAVED global pointers (bmy, 4/29/16) - INTEGER, POINTER :: Hg0_Id_List(:) => NULL() - INTEGER, POINTER :: Hg2_Id_List(:) => NULL() - INTEGER, POINTER :: HgP_Id_List(:) => NULL() - CONTAINS !EOC !------------------------------------------------------------------------- @@ -528,32 +504,6 @@ SUBROUTINE READ_HG2_PARTITIONING( Input_Opt, State_Grid, State_Met, & ENDDO !$OMP END PARALLEL DO -#ifdef BPCH_DIAG - !-------------------------------------------! - ! Store Fg and Fp in ND03 diagnostic. ! - !-------------------------------------------! - !$OMP PARALLEL DO & - !$OMP DEFAULT( SHARED ) & - !$OMP PRIVATE( I, J, L ) - DO L = 1, State_Grid%NZ - DO J = 1, State_Grid%NY - DO I = 1, State_Grid%NX - - ! error check - IF (Fg(I,J,L)+Fp(I,J,L) > 1) THEN - PRINT*, ' Fg + Fp > 1 @ ocean_mercury_mod.F90' - CALL GEOS_CHEM_STOP - ENDIF - - ENDDO - ENDDO - ENDDO - !$OMP END PARALLEL DO - - ! Increment diagnostic timestep counter. - CALL SET_Hg2_DIAG( INCREMENT=.TRUE. ) -#endif - ! Free pointers ARRAYso4 => NULL() ARRAYnit => NULL() @@ -602,7 +552,7 @@ SUBROUTINE DELIVER_SNOW_HG( SNOW_Hg2aq, I, J, State_Met, & ! ! !OUTPUT PARAMETERS: ! - REAL(fpp), INTENT(OUT) :: SNOW_Hg2aq(N_Hg_CATS) ! Hg2 from snow + REAL(fpp), INTENT(OUT) :: SNOW_Hg2aq ! Hg2 from snow ! ! !REMARKS: ! Snowpack Hg is delivered to the ocean instantaneously in an ionic pulse @@ -626,10 +576,10 @@ SUBROUTINE DELIVER_SNOW_HG( SNOW_Hg2aq, I, J, State_Met, & LOGICAL :: IS_OPEN_OCEAN, IS_OPEN_LAND ! Pointers - REAL(fpp), POINTER :: SNOW_HG_OC(:,:,:) - REAL(fpp), POINTER :: SNOW_HG_LN(:,:,:) - REAL(fpp), POINTER :: SNOW_HG_STORED_OC(:,:,:) - REAL(fpp), POINTER :: SNOW_HG_STORED_LN(:,:,:) + REAL(fpp), POINTER :: SNOW_HG_OC(:,:) + REAL(fpp), POINTER :: SNOW_HG_LN(:,:) + REAL(fpp), POINTER :: SNOW_HG_STORED_OC(:,:) + REAL(fpp), POINTER :: SNOW_HG_STORED_LN(:,:) !================================================================= ! DELIVER_SNOW_HG begins here! @@ -657,19 +607,21 @@ SUBROUTINE DELIVER_SNOW_HG( SNOW_Hg2aq, I, J, State_Met, & IS_MELT = ( State_Met%TS(I,J) >= 276e+0_fpp ) - DO NN = 1, N_Hg_CATS + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) - IS_SNOWHG_OC = ( (SNOW_HG_OC(I,J,NN) > 0e+0_fpp) .OR. & - (SNOW_HG_STORED_OC(I,J,NN) > 0e+0_fpp) ) - IS_SNOWHG_LN = ( (SNOW_HG_LN(I,J,NN) > 0e+0_fpp) .OR. & - (SNOW_HG_STORED_LN(I,J,NN) > 0e+0_fpp) ) + IS_SNOWHG_OC = ( (SNOW_HG_OC(I,J) > 0e+0_fpp) .OR. & + (SNOW_HG_STORED_OC(I,J) > 0e+0_fpp) ) + IS_SNOWHG_LN = ( (SNOW_HG_LN(I,J) > 0e+0_fpp) .OR. & + (SNOW_HG_STORED_LN(I,J) > 0e+0_fpp) ) ! OCEAN ! Check if melt conditions reached and snow in grid box IF ( IS_MELT .AND. IS_SNOWHG_OC .AND. IS_OPEN_OCEAN ) THEN ! Add all snow Hg to aqueous Hg2 reservoir - SNOW_Hg2aq(NN) = ( SNOW_HG_OC(I,J,NN) + SNOW_HG_STORED_OC(I,J,NN) ) + SNOW_Hg2aq = ( SNOW_HG_OC(I,J) + SNOW_HG_STORED_OC(I,J) ) !=========================================================== ! %%%%% HISTORY (aka netCDF diagnostics) %%%%% @@ -678,24 +630,22 @@ SUBROUTINE DELIVER_SNOW_HG( SNOW_Hg2aq, I, J, State_Met, & !=========================================================== IF ( State_Diag%Archive_EmisHg2snowToOcean ) THEN State_Diag%EmisHg2snowToOcean(I,J) = & - ( SNOW_HG_OC(I,J,NN) + SNOW_HG_STORED_OC(I,J,NN) ) / DT + ( SNOW_HG_OC(I,J) + SNOW_HG_STORED_OC(I,J) ) / DT ENDIF ! Zero reservoirs over ocean box, since we've added it ! all to the ocean - SNOW_HG_OC(I,J,NN) = 0e+0_fpp - SNOW_HG_STORED_OC(I,J,NN) = 0e+0_fpp + SNOW_HG_OC(I,J) = 0e+0_fpp + SNOW_HG_STORED_OC(I,J) = 0e+0_fpp ENDIF ! ocean box ! Zero reservoirs if there is exposed land IF ( IS_MELT .AND. IS_SNOWHG_LN .AND. IS_OPEN_LAND ) THEN - SNOW_HG_LN(I,J,NN) = 0e+0_fpp - SNOW_HG_STORED_LN(I,J,NN) = 0e+0_fpp + SNOW_HG_LN(I,J) = 0e+0_fpp + SNOW_HG_STORED_LN(I,J) = 0e+0_fpp ENDIF - ENDDO ! Hg Tracers - ! Free pointers SNOW_HG_OC => NULL() SNOW_HG_LN => NULL() @@ -749,7 +699,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! !OUTPUT PARAMETERS: ! ! Flux of Hg(0) from the ocean [kg/s] - REAL(fpp), INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY,N_Hg_CATS) + REAL(fpp), INTENT(OUT) :: FLUX(State_Grid%NX,State_Grid%NY) INTEGER, INTENT(OUT) :: RC ! Success or failure? ! ! !REMARKS: @@ -893,7 +843,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & LOGICAL :: IS_OCEAN_BOX CHARACTER(LEN=255) :: FILENAME INTEGER :: nAdvect, NA - INTEGER :: I, J, NN, C + INTEGER :: I, J, C INTEGER :: N, N_tot_oc INTEGER :: NEXTMONTH, THISMONTH INTEGER :: THISYEAR @@ -905,8 +855,8 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fpp) :: Hg2_RED, Hg2_GONE, Hg2_CONV REAL(fpp) :: FRAC_L, FRAC_O, H, TOTDEP REAL(fpp) :: oldMLD, XTAU, TOTDEPall - REAL(fpp) :: FUP(State_Grid%NX,State_Grid%NY,N_Hg_CATS) - REAL(fpp) :: FDOWN(State_Grid%NX,State_Grid%NY,N_Hg_CATS) + REAL(fpp) :: FUP(State_Grid%NX,State_Grid%NY) + REAL(fpp) :: FDOWN(State_Grid%NX,State_Grid%NY) REAL(fpp) :: X, Y, D REAL(fpp) :: NPP_tot, A_ocean, NPP_avg, RADz REAL(fpp) :: EC @@ -919,7 +869,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fpp) :: FRAC_OPEN_OCEAN, FRAC_OCEAN_OR_ICE REAL(fpp) :: OLDFLOW, RIVER_HG REAL(fpp) :: FRAC_REDUCIBLE, UVI_RATIO - REAL(fpp) :: SNOW_Hg2aq(N_Hg_CATS) + REAL(fpp) :: SNOW_Hg2aq ! Parameters REAL(fpp), PARAMETER :: EC_w = 0.0145e+0_fpp @@ -967,9 +917,9 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Pointers TYPE(SpcConc), POINTER :: Spc(:) - REAL(fpp), POINTER :: Hg0aq(:,:,:) - REAL(fpp), POINTER :: Hg2aq(:,:,:) - REAL(fpp), POINTER :: HgPaq(:,:,:) + REAL(fpp), POINTER :: Hg0aq(:,:) + REAL(fpp), POINTER :: Hg2aq(:,:) + REAL(fpp), POINTER :: HgPaq(:,:) ! Pointers to fields in the HEMCO data structure REAL(f4), POINTER :: TOMS (:,:) ! O3 @@ -1006,11 +956,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & HgPaq => State_Chm%OceanHgP ! Loop limit for use below - IF ( LSPLIT ) THEN - N_tot_oc = 2 - ELSE - N_tot_oc = 1 - ENDIF + N_tot_oc = 1 ! Molecular weight of Hg (applicable to all tagged tracers) MHg = State_Chm%SpcData(1)%Info%MW_g * 1e-3_fpp @@ -1021,19 +967,6 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Get current year to check if leap year (jaf, 8/12/11) THISYEAR = GET_YEAR() -!%% Tagged Hg simulation no longer works, so comment out (bmy, 04 Apr 2022) -!%% !----------------------------------------------- -!%% ! Check tagged & total sums (if necessary) -!%% !----------------------------------------------- -!%% IF ( USE_CHECKS .and. LSPLIT ) THEN -!%% CALL CHECK_ATMOS_MERCURY( State_Chm, State_Grid, & -!%% 'start of OCEAN_MERCURY_FLUX' ) -!%% CALL CHECK_OCEAN_MERCURY( State_Chm, State_Grid, & -!%% 'start of OCEAN_MERCURY_FLUX' ) -!%% CALL CHECK_OCEAN_FLUXES ( State_Grid, & -!%% 'start of OCEAN_MERCURY_FLUX' ) -!%% ENDIF - !----------------------------------------------------------------- ! Read O3 data from HEMCO !----------------------------------------------------------------- @@ -1147,7 +1080,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PARALLEL DO & !$OMP DEFAULT( SHARED ) & !$OMP PRIVATE( I, vi, A_M2, Hg2_RED ) & - !$OMP PRIVATE( J, NN, k_ox, OC_tot, Hg2_CONV ) & + !$OMP PRIVATE( J, k_ox, OC_tot, Hg2_CONV ) & !$OMP PRIVATE( N, TK, CHg0, k_red_bio ) & !$OMP PRIVATE( C, TC, RADz, Hg0_OX, k_red_rad ) & !$OMP PRIVATE( D, EC, k_red, OLDMLD, TOTDEPall ) & @@ -1238,16 +1171,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & State_Diag ) ! Loop over total Hg (and ocean Hg if necessary) - DO C = 1, N_tot_oc - - ! Get Hg category # - IF ( C == 1 ) NN = ID_Hg_tot - IF ( C == 2 ) NN = ID_Hg_ocn - - ! Add snowpack Hg to ocean reservoirs - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + SNOW_Hg2aq(NN) - - ENDDO + Hg2aq(I,J) = Hg2aq(I,J) + SNOW_Hg2aq !=========================================================== ! Make sure we are in an ocean box @@ -1458,41 +1382,32 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Use CDEEPATL to scale deepwater in NAtlantic IF ( UPVEL(I,J) > 0e+0_fpp ) THEN - ! Loop over total Hg (and ocean Hg if necessary) - DO C = 1, N_tot_oc - - ! Move to start of loop for use earlier (jaf, 12/8/11) - ! Grid-box latitude [degrees] - !Y = State_Grid%YMid(I,J) + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) ! Grid box longitude [degrees] X = State_Grid%XMid(I,J) - ! Get Hg category # - IF ( C == 1 ) NN = ID_Hg_tot - !-------------------------------------------------------- ! Atlantic !-------------------------------------------------------- IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= -25.0 .and. Y < 55.0 ) ) THEN !(anls,100114) - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_atl - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepatl(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepatl(2) ) ! Hg particulate !IF ( C == 1 ) THEN ! HgC(I,J) = HgC(I,J) + UPVEL(I,J) & ! * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepatl(3) ) - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepatl(3) ) !ENDIF @@ -1502,20 +1417,17 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ELSE IF ( ( X >= -180.0 .and. X < -80.0 ) .and. & ( Y >= 30.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_npa - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(2) ) ! Hg particulate !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(3) ) !ENDIF @@ -1523,20 +1435,17 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ELSE IF ( ( X >= 25.0 .and. X < 180.0 ) .and. & ( Y >= 30.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_npa - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(2) ) ! Hg particulate & !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnpa(3) ) !ENDIF @@ -1546,20 +1455,17 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ELSE IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= 55.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_nat - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnat(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnat(2) ) ! Hg particulate !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepnat(3) ) !ENDIF @@ -1569,20 +1475,17 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ELSE IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= -65.0 .and. Y < -25.0 ) ) THEN !(anls,100114) - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_sat - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepsat(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepsat(2) ) ! Hg particulate !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepsat(3) ) !ENDIF @@ -1591,20 +1494,17 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !-------------------------------------------------------- ELSE IF ( Y >= -90.0 .and. Y < -65.0 ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_ant - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepant(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepant(2) ) ! Hg particulate !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeepant(3) ) !ENDIF @@ -1613,47 +1513,36 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !-------------------------------------------------------- ELSE IF ( Y >= 70.0 .and. Y < 90.0 ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_arc - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeeparc(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeeparc(2) ) ! Hg particulate !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeeparc(3) ) !ENDIF ELSE - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_ocn - ! Hg0 (kg) - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + UPVEL(I,J) & + Hg0aq(I,J) = Hg0aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeep(1) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + UPVEL(I,J) & + Hg2aq(I,J) = Hg2aq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeep(2) ) ! Hg particulate - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + UPVEL(I,J) & + HgPaq(I,J) = HgPaq(I,J) + UPVEL(I,J) & * ( MHg * A_M2 * FRAC_O * DTSRCE * CDeep(3) ) - !ENDIF ENDIF - ENDDO - - !---------------------------------------------------------- ! Physical transport for TOTAL TRACERS, Part III: ! Downward current transport (Ekman pumping) @@ -1662,26 +1551,25 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !---------------------------------------------------------- ELSE - ! Loop over all types of tagged tracers - DO NN = 1, N_Hg_CATS + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) & + Hg0aq(I,J) = Hg0aq(I,J) & * ( 1e+0_fpp + UPVEL(I,J) * DTSRCE / & ( MLDcm * 1e-2_fpp ) ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) & + Hg2aq(I,J) = Hg2aq(I,J) & * ( 1e+0_fpp + UPVEL(I,J) * DTSRCE / & ( MLDcm * 1e-2_fpp ) ) ! Hg particulate - HgPaq(I,J,NN) = HgPaq(I,J,NN) & + HgPaq(I,J) = HgPaq(I,J) & * ( 1e+0_fpp + UPVEL(I,J) * DTSRCE / & ( MLDcm * 1e-2_fpp ) ) - ENDDO - ENDIF !=========================================================== @@ -1694,13 +1582,14 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! and NN is the Hg category # (for Hg0aq, Hg2aq, HgP) !=========================================================== - ! Loop over all Hg categories - DO NN = 1, N_Hg_CATS + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) ! Reset flux each timestep - FLUX(I,J,NN) = 0e+0_fpp - FUP(I,J,NN) = 0e+0_fpp - FDOWN(I,J,NN) = 0e+0_fpp + FLUX(I,J) = 0e+0_fpp + FUP(I,J) = 0e+0_fpp + FDOWN(I,J) = 0e+0_fpp !-------------------------------------------------------- ! Calculate new Hg(II) mass @@ -1709,7 +1598,7 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Convert to kg/box/timestep and add to Hg2aq !-------------------------------------------------------- IF ( (LArcticRiv) .AND. (Y >= 70) ) THEN - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + RIVER_HG * & + Hg2aq(I,J) = Hg2aq(I,J) + RIVER_HG * & A_M2 * DTSRCE * FRAC_O !----------------------------------------------------- @@ -1730,13 +1619,13 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Before 11/3/2009 (cdh, hamos) !! Total Hg(II) deposited on ocean surface [kg] - !TOTDEP = (WD_Hg2(I,J,NN) + DD_Hg2(I,J,NN))*FRAC_O + !TOTDEP = (WD_Hg2(I,J) + DD_Hg2(I,J))*FRAC_O ! ! Total Hg(II) deposited on ocean surface [kg] ! Includes gaseous and particulate reactive Hg(II) ! plus anthropogenic primary Hg(p) (cdh, hamos 11/3/2009) - TOTDEPall = (WD_Hg2(I,J,NN) + DD_Hg2(I,J,NN) + & - WD_HgP(I,J,NN) + DD_HgP(I,J,NN) ) + TOTDEPall = (WD_Hg2(I,J) + DD_Hg2(I,J) + & + WD_HgP(I,J) + DD_HgP(I,J) ) ! When distinguishing open ocean vs ice, only deposition ! to ocean should be included in Hg exchange (jaf, 11/28/11) @@ -1744,9 +1633,9 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & TOTDEP = TOTDEPall * FRAC_OPEN_OCEAN ! Add deposited Hg(II) to the Hg(II)tot ocean mass [kg] - Hg2aq_tot = Hg2aq(I,J,NN) + HgPaq(I,J,NN) + TOTDEP + Hg2aq_tot = Hg2aq(I,J) + HgPaq(I,J) + TOTDEP - Hg2aq(I,J,NN) = Hg2aq_tot * Frac_Hg2 + Hg2aq(I,J) = Hg2aq_tot * Frac_Hg2 ! Mass of Hg(II) --> Hg(0) !--------------------------- @@ -1816,46 +1705,46 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Now use new FRAC_REDUCIBLE for photo-reduction only. ! Retain 40% for biological reduction - !Hg2_RED = Hg2aq(I,J,NN) * 0.4d0 * k_red * DTSRCE - Hg2_RED = Hg2aq(I,J,NN) * DTSRCE * ( 0.4e+0_fpp * & + !Hg2_RED = Hg2aq(I,J) * 0.4d0 * k_red * DTSRCE + Hg2_RED = Hg2aq(I,J) * DTSRCE * ( 0.4e+0_fpp * & k_red_bio + FRAC_REDUCIBLE * k_red_rad ) ! Mass of Hg(0) --> Hg(II) - Hg0_OX = Hg0aq(I,J,NN) * k_ox * DTSRCE + Hg0_OX = Hg0aq(I,J) * k_ox * DTSRCE ! Amount of Hg(II) that is lost [kg] Hg2_GONE = Hg2_RED - Hg0_OX ! Cap Hg2_GONE with available Hg2 - IF ( Hg2_GONE > Hg2aq(I,J,NN) ) THEN - Hg2_GONE = MIN( Hg2_GONE, Hg2aq(I,J,NN) ) + IF ( Hg2_GONE > Hg2aq(I,J) ) THEN + Hg2_GONE = MIN( Hg2_GONE, Hg2aq(I,J) ) ENDIF - IF ( (Hg2_GONE * (-1e+0_fpp)) > Hg0aq(I,J,NN)) THEN - Hg2_GONE = (Hg0aq(I,J,NN)*(-1)) - !MAX (Hg2_GONE ,(Hg0aq(I,J,NN)*(-1e+0_fpp))) + IF ( (Hg2_GONE * (-1e+0_fpp)) > Hg0aq(I,J)) THEN + Hg2_GONE = (Hg0aq(I,J)*(-1)) + !MAX (Hg2_GONE ,(Hg0aq(I,J)*(-1e+0_fpp))) ENDIF ! Hg(II) ocean mass after reduction and conversion [kg] - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) - Hg2_GONE + Hg2aq(I,J) = Hg2aq(I,J) - Hg2_GONE !-------------------------------------------------------- ! Calculate new Hg(P) mass !-------------------------------------------------------- ! HgP ocean mass after conversion - HgPaq(I,J,NN) = Hg2aq_tot * ( 1 - Frac_Hg2) + HgPaq(I,J) = Hg2aq_tot * ( 1 - Frac_Hg2) !---------------------------------------------------- ! Conversion between OC and Hg !---------------------------------------------------- ! Hg/C ratio based on HgP(kg) and Stock of organic C(kg) ! HgPaq_sunk funtion of C sunk and HgP/C ratio - HgPaq_SUNK(I,J,NN) = JorgC_kg * ( HgPaq(I,J,NN ) / OC_tot_kg ) + HgPaq_SUNK(I,J) = JorgC_kg * ( HgPaq(I,J ) / OC_tot_kg ) ! HgP ocean mass after sinking [kg] - HgPaq(I,J,NN) = HgPaq(I,J,NN) - HgPaq_SUNK(I,J,NN) - HgPaq(I,J,NN) = MAX ( HgPaq(I,J,NN) , 0e+0_fpp ) + HgPaq(I,J) = HgPaq(I,J) - HgPaq_SUNK(I,J) + HgPaq(I,J) = MAX ( HgPaq(I,J) , 0e+0_fpp ) !----------------------------------------------------- ! %%%%% HISTORY DIAGNOSTIC %%%%% @@ -1870,19 +1759,16 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Calculate new Hg(0) mass !-------------------------------------------------------- - ! Hg0 tracer number (for Spc) - N = Hg0_Id_List(NN) - ! Add converted Hg(II) and subtract converted Hg(0) mass ! to the ocean mass of Hg(0) [kg] - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + Hg2_GONE + Hg0aq(I,J) = Hg0aq(I,J) + Hg2_GONE !-------------------------------------------------------- ! Calculate oceanic and gas-phase concentration of Hg(0) !-------------------------------------------------------- ! Concentration of Hg(0) in the ocean [ng/L] - CHg0aq = ( Hg0aq(I,J,NN) * 1e+11_fpp ) / & + CHg0aq = ( Hg0aq(I,J) * 1e+11_fpp ) / & ( A_M2 * FRAC_O ) / MLDcm ! Gas phase Hg(0) concentration: convert [kg] -> [ng/L] @@ -1894,34 +1780,34 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !-------------------------------------------------------- ! Compute ocean flux of Hg0 [cm/h*ng/L] - FLUX(I,J,NN) = Kw * ( CHg0aq - ( CHg0 / H ) ) + FLUX(I,J) = Kw * ( CHg0aq - ( CHg0 / H ) ) ! TURN OFF EVASION - !FLUX(I,J,NN)= MIN(0.,FLUX(I,J,NN)) + !FLUX(I,J)= MIN(0.,FLUX(I,J)) !Prior to 09 Nov 2011, H Amos --------------------- !Extra diagnostic: compute flux up and flux down - !FUP(I,J,NN) = ( Kw * CHg0aq ) - !FDOWN(I,J,NN) = ( Kw * CHg0 / H ) + !FUP(I,J) = ( Kw * CHg0aq ) + !FDOWN(I,J) = ( Kw * CHg0 / H ) !-------------------------------------------------- ! Convert [cm/h*ng/L] --> [kg/m2/s] --> [kg/s] ! Also account for ocean fraction of grid box - !FLUX(I,J,NN) = FLUX(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O + !FLUX(I,J) = FLUX(I,J) * TO_KGM2S * A_M2 * FRAC_O ! Assume fast horizontal equilibration w/in the grid box ! therefore no need to scale by open ocean fraction ! (jaf, 6/22/11) ! Evasive flux only if there is some open ocean in the ! grid box (not 100% sea ice) IF ( FRAC_OPEN_OCEAN > 0d0 ) THEN - FLUX(I,J,NN) = FLUX(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O + FLUX(I,J) = FLUX(I,J) * TO_KGM2S * A_M2 * FRAC_O ELSE - FLUX(I,J,NN) = 0e+0_fpp + FLUX(I,J) = 0e+0_fpp ENDIF !Prior to 09 Nov 2011, H Amos --------------------- - !FUP(I,J,NN) = FUP(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O - !FDOWN(I,J,NN) = FDOWN(I,J,NN) * TO_KGM2S * A_M2 * FRAC_O + !FUP(I,J) = FUP(I,J) * TO_KGM2S * A_M2 * FRAC_O + !FDOWN(I,J) = FDOWN(I,J) * TO_KGM2S * A_M2 * FRAC_O !-------------------------------------------------- !-------------------------------------------------------- @@ -1930,42 +1816,40 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !Prior to 09 Nov 2011, H Amos --------------------- ! Cap the flux w/ the available Hg(0) ocean mass - !IF ( FLUX(I,J,NN) * DTSRCE > Hg0aq(I,J,NN) ) THEN - ! FLUX(I,J,NN) = Hg0aq(I,J,NN) / DTSRCE - ! FUP(I,J,NN) = FLUX(I,J,NN)-FDOWN(I,J,NN) + !IF ( FLUX(I,J) * DTSRCE > Hg0aq(I,J) ) THEN + ! FLUX(I,J) = Hg0aq(I,J) / DTSRCE + ! FUP(I,J) = FLUX(I,J)-FDOWN(I,J) !ENDIF - IF ( FLUX(I,J,NN) * DTSRCE > Hg0aq(I,J,NN) ) THEN - FLUX(I,J,NN) = Hg0aq(I,J,NN) / DTSRCE + IF ( FLUX(I,J) * DTSRCE > Hg0aq(I,J) ) THEN + FLUX(I,J) = Hg0aq(I,J) / DTSRCE ENDIF !-------------------------------------------------- ! Cap the neg flux w/ the available Hg(0) atm mass - IF ( (-FLUX(I,J,NN) * DTSRCE ) > Spc(N)%Conc(I,J,1) ) THEN - FLUX(I,J,NN) = -Spc(N)%Conc(I,J,1) / DTSRCE + IF ( (-FLUX(I,J) * DTSRCE ) > Spc(N)%Conc(I,J,1) ) THEN + FLUX(I,J) = -Spc(N)%Conc(I,J,1) / DTSRCE ENDIF ! Cap FDOWN with available Hg(0) atm mass - !IF ((FDOWN(I,J,NN)*DTSRCE)>Spc(N)%Conc(I,J,1)) THEN - ! FDOWN(I,J,NN) = Spc(N)%Conc(I,J,1) / DTSRCE + !IF ((FDOWN(I,J)*DTSRCE)>Spc(N)%Conc(I,J,1)) THEN + ! FDOWN(I,J) = Spc(N)%Conc(I,J,1) / DTSRCE !ENDIF ! make sure Fup and Fdown do not underflow either ! debug 2x2.5 diagnostic? - FUP(I,J,NN) = MAX (FUP(I,J,NN), SMALLNUM ) - FDOWN(I,J,NN) = MAX (FDOWN(I,J,NN),SMALLNUM ) + FUP(I,J) = MAX (FUP(I,J), SMALLNUM ) + FDOWN(I,J) = MAX (FDOWN(I,J),SMALLNUM ) !-------------------------------------------------------- ! Remove amt of Hg(0) that is leaving the ocean [kg] !-------------------------------------------------------- - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) - ( FLUX(I,J,NN) * DTSRCE ) + Hg0aq(I,J) = Hg0aq(I,J) - ( FLUX(I,J) * DTSRCE ) ! Make sure Hg0aq does not underflow (cdh, bmy, 3/28/06) - Hg0aq(I,J,NN) = MAX( Hg0aq(I,J,NN), SMALLNUM ) + Hg0aq(I,J) = MAX( Hg0aq(I,J), SMALLNUM ) - !Hgaq_tot = HgC(I,J) + Hg0aq(I,J,NN) + Hg2aq(I,J,NN) - Hgaq_tot(I,J,NN) = HgPaq(I,J,NN) + Hg0aq(I,J,NN) + Hg2aq(I,J,NN) - - ENDDO + !Hgaq_tot = HgC(I,J) + Hg0aq(I,J) + Hg2aq(I,J) + Hgaq_tot(I,J) = HgPaq(I,J) + Hg0aq(I,J) + Hg2aq(I,J) !----------------------------------------------------------- ! %%%%% HISTORY (aka netCDF) DIAGNOSTICS %%%%% @@ -1975,48 +1859,45 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & ! NOTE: We have converted units from kg !----------------------------------------------------------- - ! NOTE: for now just use the total category - NN = 1 - ! Mass of oceanic Hg0 [kg] IF ( State_Diag%Archive_MassHg0inOcean ) THEN - State_Diag%MassHg0inOcean(I,J) = Hg0aq(I,J,NN) + State_Diag%MassHg0inOcean(I,J) = Hg0aq(I,J) ENDIF ! Mass of oceanic Hg2 [kg] IF ( State_Diag%Archive_MassHg2inOcean ) THEN - State_Diag%MassHg2inOcean(I,J) = Hg2aq(I,J,NN) + State_Diag%MassHg2inOcean(I,J) = Hg2aq(I,J) ENDIF ! Mass of oceanic HgP [kg] IF ( State_Diag%Archive_MassHgPinOcean ) THEN - State_Diag%MassHgPinOcean(I,J) = HgPaq(I,J,NN) + State_Diag%MassHgPinOcean(I,J) = HgPaq(I,J) ENDIF ! Total oceanic mercury [kg] IF ( State_Diag%Archive_MassHgTotalInOcean ) THEN - State_Diag%MassHgTotalInOcean(I,J) = Hgaq_tot(I,J,NN) + State_Diag%MassHgTotalInOcean(I,J) = Hgaq_tot(I,J) ENDIF ! Flux of Hg2 sunk to deep ocean [kg/s] IF ( State_Diag%Archive_FluxHg2toDeepOcean ) THEN - State_Diag%FluxHg2toDeepOcean(I,J) = HgPaq_SUNK(I,J,NN) / DTSRCE + State_Diag%FluxHg2toDeepOcean(I,J) = HgPaq_SUNK(I,J) / DTSRCE ENDIF ! Ocean-to-air and air-to-ocean fluxes - IF ( FLUX(I,J,NN) > 0e+0_fpp ) THEN + IF ( FLUX(I,J) > 0e+0_fpp ) THEN ! Volatization flux of Hg0 from ocean to air ! NOTE: Units are kg/s IF ( State_Diag%Archive_FluxHg0fromOceantoAir ) THEN - State_Diag%FluxHg0fromOceanToAir(I,J) = FLUX(I,J,NN) + State_Diag%FluxHg0fromOceanToAir(I,J) = FLUX(I,J) ENDIF - ELSE IF ( FLUX(I,J,NN) < 0e+0_fpp) THEN + ELSE IF ( FLUX(I,J) < 0e+0_fpp) THEN ! Drydep flux of Hg0 from air to ocean IF ( State_Diag%Archive_FluxHg0fromAirToOcean ) THEN - State_Diag%FluxHg0fromAirToOcean(I,J) = ABS( FLUX(I,J,NN) ) + State_Diag%FluxHg0fromAirToOcean(I,J) = ABS( FLUX(I,J) ) ENDIF ENDIF @@ -2032,11 +1913,12 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & !============================================================== ELSE - DO NN = 1, N_Hg_CATS - FLUX(I,J,NN) = 0e+0_fpp - FUP(I,J,NN)=0e+0_fpp - FDOWN(I,J,NN)=0e+0_fpp - ENDDO + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) + FLUX(I,J) = 0e+0_fpp + FUP(I,J)=0e+0_fpp + FDOWN(I,J)=0e+0_fpp ENDIF @@ -2057,21 +1939,6 @@ SUBROUTINE OCEAN_MERCURY_FLUX( Input_Opt, State_Chm, State_Diag, State_Grid, & TOMS_PD => NULL() TOMS_LT => NULL() -!%% Tagged simulation no longer works so comment out (bmy, 04 Apr 2022) -!%% !================================================================= -!%% ! Check tagged & total sums (if necessary) -!%% !================================================================= -!%% IF ( USE_CHECKS .and. LSPLIT ) THEN -!%% CALL CHECK_ATMOS_MERCURY( State_Chm, State_Grid, & -!%% 'end of OCEAN_MERCURY_FLUX' ) -!%% CALL CHECK_OCEAN_MERCURY( State_Chm, State_Grid, & -!%% 'end of OCEAN_MERCURY_FLUX' ) -!%% CALL CHECK_OCEAN_FLUXES ( State_Grid, & -!%% 'end of OCEAN_MERCURY_FLUX' ) -!%% CALL CHECK_FLUX_OUT( State_Grid, FLUX, & -!%% 'end of OCEAN_MERCURY_FLUX' ) -!%% ENDIF - END SUBROUTINE OCEAN_MERCURY_FLUX !EOC !------------------------------------------------------------------------------ @@ -2390,7 +2257,7 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ! ! !LOCAL VARIABLES: ! - INTEGER :: C, NN, N_tot_oc + INTEGER :: C, NN INTEGER :: K, L REAL(fpp) :: A_M2, DELTAH, FRAC_O, MHg REAL(fpp) :: X, Y !(added anls 01/05/09) @@ -2399,9 +2266,9 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & LOGICAL :: LSPLIT ! Pointers - REAL(fpp), POINTER :: Hg0aq(:,:,:) - REAL(fpp), POINTER :: Hg2aq(:,:,:) - REAL(fpp), POINTER :: HgPaq(:,:,:) + REAL(fpp), POINTER :: Hg0aq(:,:) + REAL(fpp), POINTER :: Hg2aq(:,:) + REAL(fpp), POINTER :: HgPaq(:,:) !================================================================= ! MLD_ADJUSTMENT begins here! @@ -2415,13 +2282,6 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & Hg2aq => State_Chm%OceanHg2 HgPaq => State_Chm%OceanHgP - ! Loop limit for use below - IF ( LSPLIT ) THEN - N_tot_oc = 2 - ELSE - N_tot_oc = 1 - ENDIF - ! Grid box surface area [m2] A_M2 = State_Grid%Area_M2(I,J) @@ -2435,7 +2295,7 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & FRAC_O = State_Met%FROCEAN(I,J) ! Molecular weight of Hg (valid for all tagged tracers) - MHg = State_Chm%SpcData(ID_Hg_tot)%Info%MW_g * 1e-3_fpp + MHg = State_Chm%SpcData(1)%Info%MW_g * 1e-3_fpp ! Test if MLD increased IF ( MLDnew > MLDold ) THEN @@ -2461,11 +2321,9 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ! Grid box longitude [degrees] X = State_Grid%XMid(I,J) - ! Loop over total Hg (and ocean Hg if necessary) - DO C = 1, N_tot_oc - - ! Get Hg category number - IF ( C == 1 ) NN = ID_Hg_tot + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) !------------------------------------------------ ! Atlantic @@ -2473,20 +2331,17 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= -25.0 .and. Y < 55.0 ) ) THEN !(anls,100114) - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_atl - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepatl(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepatl(2) * MHg * A_M2 * FRAC_O ) ! HgP !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepatl(3) * MHg * A_M2 * FRAC_O ) !ENDIF @@ -2496,20 +2351,17 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ELSE IF ( ( X >= -180.0 .and. X < -80.0 ) .and. & ( Y >= 30.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_npa - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepnpa(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepnpa(2) * MHg * A_M2 * FRAC_O ) ! HgP !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepnpa(3) * MHg * A_M2 * FRAC_O ) !ENDIF @@ -2519,22 +2371,17 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ELSE IF ( ( X >= 25.0 .and. X < 180.0 ) .and. & ( Y >= 30.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_npa - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepnpa(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepnpa(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepnpa(3) * MHg * A_M2 * FRAC_O ) - !ENDIF !------------------------------------------------ ! North Atlantic @@ -2542,20 +2389,16 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ELSE IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= 55.0 .and. Y < 70.0 ) ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_nat - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepnat(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepnat(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepnat(3) * MHg * A_M2 * FRAC_O ) !ENDIF @@ -2565,93 +2408,72 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ELSE IF ( ( X >= -80.0 .and. X < 25.0 ) .and. & ( Y >= -65.0 .and. Y < -25.0 ) ) THEN !(anls,100114) - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_sat - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepsat(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepsat(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepsat(3) * MHg * A_M2 * FRAC_O ) - !ENDIF !------------------------------------------------ ! Antarctic !------------------------------------------------ ELSE IF ( Y >= -90.0 .and. Y < -65.0 ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_ant - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeepant(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeepant(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeepant(3) * MHg * A_M2 * FRAC_O ) - !ENDIF !------------------------------------------------ ! Arctic !------------------------------------------------ ELSE IF ( Y >= 70.0 .and. Y < 90.0 ) THEN - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_arc - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeeparc(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeeparc(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeeparc(3) * MHg * A_M2 * FRAC_O ) - !ENDIF ELSE - ! Assign region tag - IF ( C == 2 ) NN = ID_Hg_ocn - ! Hg0 - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) + & + Hg0aq(I,J) = Hg0aq(I,J) + & ( DELTAH * CDeep(1) * MHg * A_M2 * FRAC_O ) ! Hg2 - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) + & + Hg2aq(I,J) = Hg2aq(I,J) + & ( DELTAH * CDeep(2) * MHg * A_M2 * FRAC_O ) ! HgP - !IF ( C == 1 ) THEN - HgPaq(I,J,NN) = HgPaq(I,J,NN) + & + HgPaq(I,J) = HgPaq(I,J) + & ( DELTAH * CDeep(3) * MHg * A_M2 * FRAC_O ) - !ENDIF ENDIF - ENDDO ELSE !============================================================== - ! IF MIXED LAYER DEPTH HAS DECREASED: + ! IF MIXED LAYER DEPTH HAS CREASED: ! ! Conserve concentration, but shed mass for ALL tracers. ! Mass changes by same ratio as volume. @@ -2660,12 +2482,12 @@ SUBROUTINE MLD_ADJUSTMENT( I, J, MLDold, MLDnew, Input_Opt, & ! Avoid dividing by zero IF ( MLDold > 0e+0_fpp ) THEN - ! Update Hg0 and Hg2 categories - DO NN = 1, N_Hg_CATS - Hg0aq(I,J,NN) = Hg0aq(I,J,NN) * ( MLDnew / MLDold ) - Hg2aq(I,J,NN) = Hg2aq(I,J,NN) * ( MLDnew / MLDold ) - HgPaq(I,J,NN) = HgPaq(I,J,NN) * ( MLDnew / MLDold ) - ENDDO + !%%% Removed loop over Hg categories here + !%%% and also removed NN as the 3rd array index + !%%% -- Bob Yantosca (23 Jun 2022) + Hg0aq(I,J) = Hg0aq(I,J) * ( MLDnew / MLDold ) + Hg2aq(I,J) = Hg2aq(I,J) * ( MLDnew / MLDold ) + HgPaq(I,J) = HgPaq(I,J) * ( MLDnew / MLDold ) ENDIF @@ -2683,543 +2505,6 @@ END SUBROUTINE MLD_ADJUSTMENT !------------------------------------------------------------------------------ !BOP ! -! !IROUTINE: check_atmos_mercury -! -! !DESCRIPTION: Subroutine CHECK\_ATMOS\_MERCURY tests whether the total and -! tagged tracers the GEOS-CHEM species array sum properly within each grid box. -! (cdh, bmy, 3/28/06) -!\\ -!\\ -! !INTERFACE: -! - SUBROUTINE CHECK_ATMOS_MERCURY( State_Chm, State_Grid, LOC ) -! -! !USES: -! - USE ERROR_MOD, ONLY : ERROR_STOP - USE Species_Mod, ONLY : SpcConc - USE State_Chm_Mod, ONLY : ChmState - USE State_Grid_Mod, ONLY : GrdState -! -! !INPUT PARAMETERS: -! - CHARACTER(LEN=*), INTENT(IN) :: LOC ! Calling routine - TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object - TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object -! -! !REVISION HISTORY: -! See https://github.com/geoschem/geos-chem for complete history -!EOP -!------------------------------------------------------------------------------ -!BOC -! -! !LOCAL VARIABLES: -! - LOGICAL :: FLAG - INTEGER :: I, J, L - INTEGER :: N, NN - REAL(fpp) :: Hg0_tot, Hg0_tag, RELERR0, ABSERR0 - REAL(fpp) :: Hg2_tot, Hg2_tag, RELERR2, ABSERR2 - REAL(fpp) :: HgP_tot, HgP_tag, RELERRP, ABSERRP - - ! Pointers - TYPE(SpcConc), POINTER :: Spc(:) - - !================================================================= - ! CHECK_ATMOS_MERCURY begins here! - !================================================================= - - ! Set error flags - FLAG = .FALSE. - - ! Point to chemical species array [kg] - Spc => State_Chm%Species - - ! Loop over grid boxes - !OMP PARALLEL DO & - !OMP DEFAULT( SHARED ) & - !OMP PRIVATE( I, J, L, N, NN ) & - !OMP PRIVATE( Hg0_tot, RELERR0, ABSERR0 ) & - !OMP PRIVATE( Hg2_tot, RELERR2, ABSERR2 ) & - !OMP PRIVATE( HgP_tot, RELERRP, ABSERRP ) & - !OMP REDUCTION( +: Hg0_tag, Hg2_tag, HgP_tag ) - DO L = 1, State_Grid%NZ - DO J = 1, State_Grid%NY - DO I = 1, State_Grid%NX - - ! Initialize - Hg0_tot = 0e+0_fpp - Hg0_tag = 0e+0_fpp - RELERR0 = 0e+0_fpp - ABSERR0 = 0e+0_fpp - Hg2_tot = 0e+0_fpp - Hg2_tag = 0e+0_fpp - RELERR2 = 0e+0_fpp - ABSERR2 = 0e+0_fpp - HgP_tot = 0e+0_fpp - Hgp_tag = 0e+0_fpp - RELERRP = 0e+0_fpp - ABSERRP = 0e+0_fpp - - !-------- - ! Hg(0) - !-------- - - ! Total Hg(0) - N = Hg0_Id_List(ID_Hg_tot) - Hg0_tot = Spc(N)%Conc(I,J,L) - - ! Sum of tagged Hg(0) - DO NN = 2, N_Hg_CATS - N = Hg0_Id_List(NN) - Hg0_tag = Hg0_tag + Spc(N)%Conc(I,J,L) - ENDDO - - ! Absolute error for Hg0 - ABSERR0 = ABS( Hg0_tot - Hg0_tag ) - - ! Relative error for Hg0 (avoid div by zero) - IF ( Hg0_tot > 0e+0_fpp ) THEN - RELERR0 = ABS( ( Hg0_tot - Hg0_tag ) / Hg0_tot ) - ELSE - RELERR0 = -999e+0_fpp - ENDIF - - !-------- - ! Hg(II) - !-------- - - ! Total Hg(II) - N = Hg2_Id_List(ID_Hg_tot) - Hg2_tot = Spc(N)%Conc(I,J,L) - - ! Sum of tagged Hg(II) - DO NN = 2, N_Hg_CATS - N = Hg2_Id_List(NN) - Hg2_tag = Hg2_tag + Spc(N)%Conc(I,J,L) - ENDDO - - ! Absolute error for Hg2 - ABSERR2 = ABS( Hg2_tot - Hg2_tag ) - - ! Relative error for Hg2 (avoid div by zero) - IF ( Hg2_tot > 0e+0_fpp ) THEN - RELERR2 = ABS( ( Hg2_tot - Hg2_tag ) / Hg2_tot ) - ELSE - RELERR2 = -999e+0_fpp - ENDIF - - !-------- - ! HgP - !-------- - - ! Total Hg(P) - N = HgP_Id_List(ID_Hg_tot) - HgP_tot = Spc(N)%Conc(I,J,L) - - ! Sum of tagged Hg(P) - DO NN = 2, N_Hg_CATS - N = HgP_Id_List(NN) - IF ( N > 0 ) HgP_tag = HgP_tag + Spc(N)%Conc(I,J,L) - ENDDO - - ! Absolute error for HgP - ABSERRP = ABS( HgP_tot - HgP_tag ) - - ! Relative error for HgP (avoid div by zero) - IF ( HgP_tot > 0e+0_fpp ) THEN - RELERRP = ABS( ( HgP_tot - HgP_tag ) / HgP_tot ) - ELSE - RELERRP = -999e+0_fpp - ENDIF - - !---------------------------- - ! Hg(0) error is too large - !---------------------------- - IF ( RELERR0 > MAX_RELERR .and. ABSERR0 > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 100 ) I, J, L, Hg0_tot, Hg0_tag, RELERR0, ABSERR0 - !OMP END CRITICAL - ENDIF - - !---------------------------- - ! Hg(0) error is too large - !---------------------------- - IF ( RELERR2 > MAX_RELERR .and. ABSERR2 > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 110 ) I, J, L, Hg2_tot, Hg2_tag, RELERR2, ABSERR2 - !OMP END CRITICAL - ENDIF - - !---------------------------- - ! HgP error is too large - !---------------------------- - IF ( RELERRP > MAX_RELERR .and. ABSERRP > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 120 ) I, J, L, HgP_tot, HgP_tag, RELERRP, ABSERRP - !OMP END CRITICAL - ENDIF - ENDDO - ENDDO - ENDDO - !OMP END PARALLEL DO - - ! Free pointer - Spc => NULL() - - ! FORMAT strings -100 FORMAT( 'Hg0 error: ', 3i5, 4es13.6 ) -110 FORMAT( 'Hg2 error: ', 3i5, 4es13.6 ) -120 FORMAT( 'HgP error: ', 3i5, 4es13.6 ) - - ! Stop if Hg0 and Hg2 errors are too large - IF ( FLAG ) THEN - CALL ERROR_STOP( 'Tagged Hg0, Hg2, HgP do not add up!', LOC ) - ENDIF - - END SUBROUTINE CHECK_ATMOS_MERCURY -!EOC -!------------------------------------------------------------------------------ -! GEOS-Chem Global Chemical Transport Model ! -!------------------------------------------------------------------------------ -!BOP -! -! !IROUTINE: check_ocean_mercury -! -! !DESCRIPTION: Subroutine CHECK\_OCEAN_MERCURY tests whether tagged tracers in -! Hg0aq and Hg2aq add properly within each grid box. (cdh, bmy, 3/28/06) -!\\ -!\\ -! !INTERFACE: -! - SUBROUTINE CHECK_OCEAN_MERCURY( State_Chm, State_Grid, LOC ) -! -! !USES: -! - USE ERROR_MOD, ONLY : ERROR_STOP - USE State_Chm_Mod, ONLY : ChmState - USE State_Grid_Mod, ONLY : GrdState -! -! !INPUT PARAMETERS: -! - CHARACTER(LEN=*), INTENT(IN) :: LOC ! Calling routine - TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object - TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object -! -! !REVISION HISTORY: -! See https://github.com/geoschem/geos-chem for complete history -!EOP -!------------------------------------------------------------------------------ -!BOC -! -! !LOCAL VARIABLES: -! - LOGICAL, SAVE :: FIRST = .TRUE. - LOGICAL :: FLAG - INTEGER :: I, J - REAL(fpp) :: Hg0_tot, Hg0_tag, RELERR0, ABSERR0 - REAL(fpp) :: Hg2_tot, Hg2_tag, RELERR2, ABSERR2 - - !================================================================= - ! CHECK_OCEAN_MERCURY begins here! - !================================================================= - - ! Set error condition flag - FLAG = .FALSE. - - ! Loop over ocean surface boxes - !OMP PARALLEL DO & - !OMP DEFAULT( SHARED ) & - !OMP PRIVATE( I, J ) & - !OMP PRIVATE( Hg0_tot, Hg0_tag, RELERR0, ABSERR0 ) - !OMP PRIVATE( Hg2_tot, Hg2_tag, RELERR2, ABSERR2 ) - DO J = 1, State_Grid%NY - DO I = 1, State_Grid%NX - - !-------------------------------------- - ! Relative and absolute errors for Hg0 - !-------------------------------------- - Hg0_tot = State_Chm%OceanHg0(I,J,ID_Hg_tot) - Hg0_tag = SUM( State_Chm%OceanHg0(I,J,2:N_Hg_CATS) ) - ABSERR0 = ABS( Hg0_tot - Hg0_tag ) - - ! Avoid div by zero - IF ( Hg0_tot > 0e+0_fpp ) THEN - RELERR0 = ABS( ( Hg0_tot - Hg0_tag ) / Hg0_tot ) - ELSE - RELERR0 = -999e+0_fpp - ENDIF - - !-------------------------------------- - ! Relative and absolute errors for Hg2 - !-------------------------------------- - Hg2_tot = State_Chm%OceanHg2(I,J,ID_Hg_tot) - Hg2_tag = SUM( State_Chm%OceanHg2(I,J,2:N_Hg_CATS) ) - ABSERR2 = ABS( Hg2_tot - Hg2_tag ) - - ! Avoid div by zero - IF ( Hg2_tot > 0e+0_fpp ) THEN - RELERR2 = ABS( ( Hg2_tot - Hg2_tag ) / Hg2_tot ) - ELSE - RELERR2 = -999e+0_fpp - ENDIF - - !-------------------------------------- - ! Hg(0) error is too large - !-------------------------------------- - IF ( RELERR0 > MAX_RELERR .and. ABSERR0 > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 100 ) I, J, Hg0_tot, Hg0_tag, RELERR0, ABSERR0 - !OMP END CRITICAL - ENDIF - - !-------------------------------------- - ! Hg(II) error is too large - !-------------------------------------- - IF ( RELERR2 > MAX_RELERR .and. ABSERR2 > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 110 ) I, J, Hg2_tot, Hg2_tag, RELERR2, ABSERR2 - !OMP END CRITICAL - ENDIF - - ENDDO - ENDDO - !OMP END PARALLEL DO - - ! FORMAT strings -100 FORMAT( 'Hg0aq error: ', 2i5, 4es13.6 ) -110 FORMAT( 'Hg2aq error: ', 2i5, 4es13.6 ) - - ! Stop if Hg0 and Hg2 errors are too large - IF ( FLAG ) THEN - CALL ERROR_STOP( 'Tagged Hg0aq, Hg2aq do not add up!', LOC ) - ENDIF - - END SUBROUTINE CHECK_OCEAN_MERCURY -!EOC -!------------------------------------------------------------------------------ -! GEOS-Chem Global Chemical Transport Model ! -!------------------------------------------------------------------------------ -!BOP -! -! !IROUTINE: check_ocean_fluxes -! -! !DESCRIPTION: Subroutine CHECK\_OCEAN\_FLUXES tests whether the drydep and -! wetdep fluxes in DD_Hg2 and WD_Hg2 sum together in each grid box. (cdh, bmy, -! 3/28/06) -!\\ -!\\ -! !INTERFACE: -! - SUBROUTINE CHECK_OCEAN_FLUXES( State_Grid, LOC ) -! -! !USES: -! - USE DEPO_MERCURY_MOD, ONLY : DD_Hg2, WD_Hg2, DD_HgP, WD_HgP - USE ERROR_MOD, ONLY : ERROR_STOP - USE State_Grid_Mod, ONLY : GrdState -! -! !INPUT PARAMETERS: -! - CHARACTER(LEN=*), INTENT(IN) :: LOC ! Calling routine - TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object -! -! !REVISION HISTORY: -! See https://github.com/geoschem/geos-chem for complete history -!EOP -!------------------------------------------------------------------------------ -!BOC -! -! !LOCAL VARIABLES: -! - LOGICAL :: FLAG - INTEGER :: I, J - REAL(fpp) :: DD_tot, DD_tag - REAL(fpp) :: DD_RELERR, DD_ABSERR - REAL(fpp) :: WD_tot, WD_tag - REAL(fpp) :: WD_RELERR, WD_ABSERR - - !================================================================= - ! CHECK_OCEAN_MERCURY begins here! - !================================================================= - - ! Echo - WRITE( 6, 100 ) -100 FORMAT( ' - In CHECK_OCEAN_FLUXES' ) - - ! Set error condition flag - FLAG = .FALSE. - - ! Loop over ocean surface boxes - !OMP PARALLEL DO & - !OMP DEFAULT( SHARED ) & - !OMP PRIVATE( I, J ) & - !OMP PRIVATE( DD_tot, DD_tag, DD_RELERR, DD_ABSERR ) & - !OMP PRIVATE( WD_tot, WD_tag, WD_RELERR, WD_ABSERR ) - DO J = 1, State_Grid%NY - DO I = 1, State_Grid%NX - - !--------------------------------------- - ! Absolute & relative errors for DD_Hg2 - !--------------------------------------- - DD_tot = DD_Hg2(I,J,1) - DD_tag = SUM( DD_Hg2(I,J,2:N_Hg_CATS) ) - DD_ABSERR = ABS( DD_tot - DD_tag ) - - ! Avoid div by zero - IF ( DD_tot > 0e+0_fpp ) THEN - DD_RELERR = ABS( ( DD_tot - DD_tag ) / DD_tot ) - ELSE - DD_RELERR = -999e+0_fpp - ENDIF - - !--------------------------------------- - ! Absolute & relative errors for WD_Hg2 - !--------------------------------------- - WD_tot = WD_Hg2(I,J,1) - WD_tag = SUM( WD_Hg2(I,J,2:N_Hg_CATS) ) - WD_ABSERR = ABS( WD_tot - WD_tag ) - - ! Avoid div by zero - IF ( WD_tot > 0e+0_fpp ) THEN - WD_RELERR = ABS( ( WD_tot - WD_tag ) / WD_tot ) - ELSE - WD_RELERR = -999e+0_fpp - ENDIF - - !--------------------------------------- - ! DD flux error is too large - !--------------------------------------- - IF ( DD_RELERR > MAX_RELERR .and. DD_ABSERR > MAX_FLXERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 110 ) I, J, DD_tot, DD_tag, DD_RELERR, DD_ABSERR - !OMP END CRITICAL - ENDIF - - !--------------------------------------- - ! WD flux error is too large - !--------------------------------------- - IF ( WD_RELERR > MAX_RELERR .and. WD_ABSERR > MAX_FLXERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 120 ) I, J, WD_tot, WD_tag, WD_RELERR, WD_ABSERR - !OMP END CRITICAL - ENDIF - - ENDDO - ENDDO - !OMP END PARALLEL DO - - ! FORMAT strings -110 FORMAT( 'DD_Hg2 flux error: ', 2i5, 4es13.6 ) -120 FORMAT( 'WD_Hg2 flux error: ', 2i5, 4es13.6 ) - - ! Stop if Hg0 and Hg2 errors are too large - IF ( FLAG ) THEN - CALL ERROR_STOP( 'Tagged DD, WD fluxes do not add up!', LOC ) - ENDIf - - END SUBROUTINE CHECK_OCEAN_FLUXES -!EOC -!------------------------------------------------------------------------------ -! GEOS-Chem Global Chemical Transport Model ! -!------------------------------------------------------------------------------ -!BOP -! -! !IROUTINE: check_flux_out -! -! !DESCRIPTION: Subroutine CHECK\_FLUX\_OUT tests whether tagged quantities in -! FLUX sum together in each grid box. (cdh, bmy, 3/20/06) -!\\ -!\\ -! !INTERFACE: -! - SUBROUTINE CHECK_FLUX_OUT( State_Grid, FLUX, LOC ) -! -! !USES: -! - USE ERROR_MOD, ONLY : ERROR_STOP - USE State_Grid_Mod, ONLY : GrdState -! -! !INPUT PARAMETERS: -! - TYPE(GrdState), INTENT(IN) :: State_Grid - REAL(fpp), INTENT(IN) :: FLUX(State_Grid%NX,State_Grid%NY,N_Hg_CATS) - CHARACTER(LEN=*), INTENT(IN) :: LOC - - ! Local variables - LOGICAL :: FLAG - INTEGER :: I, J - REAL(fpp) :: FLX_tot, FLX_tag - REAL(fpp) :: FLX_RELERR, FLX_ABSERR - - !================================================================= - ! CHECK_FLUX_OUT begins here! - !================================================================= - - ! Echo - WRITE( 6, 100 ) -100 FORMAT( ' - In CHECK_FLUX_OUT' ) - - ! Set error condition flag - FLAG = .FALSE. - - ! Loop over ocean surface boxes - !OMP PARALLEL DO & - !OMP DEFAULT( SHARED ) & - !OMP PRIVATE( I, J, FLX_tot, FLX_tag, FLX_err ) - DO J = 1, State_Grid%NY - DO I = 1, State_Grid%NX - - !---------------------------------------- - ! Absolute & relative errors for FLX_Hg2 - !---------------------------------------- - FLX_tot = FLUX(I,J,1) - FLX_tag = SUM( FLUX(I,J,2:N_Hg_CATS) ) - FLX_ABSERR = ABS( FLX_tot - FLX_tag ) - - ! Avoid div by zero - IF ( FLX_tot > 0e+0_fpp ) THEN - FLX_RELERR = ABS( ( FLX_tot - FLX_tag ) / FLX_tot ) - ELSE - FLX_RELERR = -999e+0_fpp - ENDIF - - !---------------------------- - ! Flux error is too large - !---------------------------- - IF ( FLX_RELERR > MAX_RELERR .and. FLX_ABSERR > MAX_ABSERR ) THEN - !OMP CRITICAL - FLAG = .TRUE. - WRITE( 6, 110 ) I, J, FLX_tot, FLX_tag, FLX_RELERR, FLX_ABSERR - !OMP END CRITICAL - ENDIF - - ENDDO - ENDDO - !OMP END PARALLEL DO - - ! FORMAT strings -110 FORMAT( 'FLX_Hg2 flux error: ', 2i5, 4es13.6 ) - - ! Stop if Hg0 and Hg2 errors are too large - IF ( FLAG ) THEN - CALL ERROR_STOP( 'Tagged emission fluxes do not add up!', LOC ) - ENDIf - - END SUBROUTINE CHECK_FLUX_OUT -!EOC -!------------------------------------------------------------------------------ -! GEOS-Chem Global Chemical Transport Model ! -!------------------------------------------------------------------------------ -!BOP -! ! !IROUTINE: init_land_mercury ! ! !DESCRIPTION: Subroutine INIT\_OCEAN\_MERCURY allocates and zeroes all @@ -3282,47 +2567,6 @@ SUBROUTINE INIT_OCEAN_MERCURY( Input_Opt, State_Chm, State_Grid, RC ) ! Exit if this is a dry-run IF ( Input_Opt%DryRun ) RETURN - ! Store the # of tagHg categories in a module variable - 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 - - ! Loop over all Hg species - DO N = 1, State_Chm%nSpecies - - ! Point to the Species Database entry for tracer # N - SpcInfo => State_Chm%SpcData(N)%Info - - ! Save ocean Hg categories into module variable indices - SELECT CASE( TRIM( SpcInfo%Name ) ) - CASE( 'Hg0' ) - ID_Hg_tot = SpcInfo%Hg_Cat - CASE( 'Hg0_atl' ) - ID_Hg_atl = SpcInfo%Hg_Cat - CASE( 'Hg0_nat' ) - ID_Hg_nat = SpcInfo%Hg_Cat - CASE( 'Hg0_sat' ) - ID_Hg_sat = SpcInfo%Hg_Cat - CASE( 'Hg0_npa' ) - ID_Hg_npa = SpcInfo%Hg_Cat - CASE( 'Hg0_arc' ) - ID_Hg_arc = SpcInfo%Hg_Cat - CASE( 'Hg0_ant' ) - ID_Hg_ant = SpcInfo%Hg_Cat - CASE( 'Hg0_ocn' ) - ID_Hg_ocn = SpcInfo%Hg_Cat - CASE DEFAULT - ! Do nothing - END SELECT - - ! Free pointer - SpcInfo => NULL() - - ENDDO - !----------------------------------------------------------------- ! Allocate these arrays regardless of whether you are using ! the dynamic ocean. These are needed for Hg2 partitioning. @@ -3385,8 +2629,7 @@ SUBROUTINE INIT_OCEAN_MERCURY( Input_Opt, State_Chm, State_Grid, RC ) Fp = 0e+0_fpp !eds 5/15/12 fix - ALLOCATE( Hgaq_tot (State_Grid%NX, State_Grid%NY, N_Hg_CATS ), & - STAT=RC ) + ALLOCATE( Hgaq_tot (State_Grid%NX, State_Grid%NY ), STAT=RC ) IF ( RC /= 0 ) CALL ALLOC_ERR( 'Hgaq_tot' ) Hgaq_tot = 0e+0_fpp @@ -3394,9 +2637,6 @@ SUBROUTINE INIT_OCEAN_MERCURY( Input_Opt, State_Chm, State_Grid, RC ) ! These variables are only needed for the dynamic ocean option !----------------------------------------------------------------- - ! Turn on error checks for tagged & total sums? - USE_CHECKS = Input_Opt%USE_CHECKS - ! Get current year for scale factor application THISYEAR = GET_YEAR() @@ -3510,8 +2750,7 @@ SUBROUTINE INIT_OCEAN_MERCURY( Input_Opt, State_Chm, State_Grid, RC ) IF ( RC /= 0 ) CALL ALLOC_ERR( 'dMLD' ) dMLD = 0e+0_fpp - ALLOCATE( HgPaq_SUNK (State_Grid%NX, State_Grid%NY, N_Hg_CATS ), & - STAT=RC ) + ALLOCATE( HgPaq_SUNK( State_Grid%NX, State_Grid%NY ), STAT=RC ) IF ( RC /= 0 ) CALL ALLOC_ERR( 'HgPaq_SUNK' ) HgPaq_SUNK = 0e+0_fpp @@ -3577,9 +2816,6 @@ SUBROUTINE CLEANUP_OCEAN_MERCURY UPVEL => NULL() dMLD1 => NULL() dMLD2 => NULL() - Hg0_Id_List => NULL() - Hg2_Id_List => NULL() - HgP_Id_List => NULL() END SUBROUTINE CLEANUP_OCEAN_MERCURY !EOC diff --git a/GeosCore/vdiff_mod.F90 b/GeosCore/vdiff_mod.F90 index 438233695..e130bd99b 100644 --- a/GeosCore/vdiff_mod.F90 +++ b/GeosCore/vdiff_mod.F90 @@ -1932,8 +1932,6 @@ SUBROUTINE VDIFFDR( Input_Opt, State_Chm, State_Diag, & USE State_Grid_Mod, ONLY : GrdState USE State_Met_Mod, ONLY : MetState USE TIME_MOD, ONLY : GET_TS_CONV, GET_TS_EMIS, GET_TS_CHEM - USE DEPO_MERCURY_MOD, ONLY : ADD_Hg2_DD, ADD_HgP_DD - USE DEPO_MERCURY_MOD, ONLY : ADD_Hg2_SNOWPACK USE OCEAN_MERCURY_MOD, ONLY : Fg !hma USE OCEAN_MERCURY_MOD, ONLY : OMMFP => Fp USE OCEAN_MERCURY_MOD, ONLY : LHg2HalfAerosol !cdh diff --git a/GeosCore/wetscav_mod.F90 b/GeosCore/wetscav_mod.F90 index 9dcb6c7fb..cda7f76f2 100644 --- a/GeosCore/wetscav_mod.F90 +++ b/GeosCore/wetscav_mod.F90 @@ -3751,23 +3751,17 @@ SUBROUTINE WETDEP( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Check if it is a gaseous Hg2 tag IF ( SpcInfo%Is_Hg2 ) THEN - ! Get the category # for this Hg2 species - Hg_Cat = SpcInfo%Hg_Cat - ! Archive wet-deposited Hg2 - CALL ADD_Hg2_WD ( I, J, Hg_Cat, DEP_HG ) - CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, DEP_HG, & + CALL ADD_Hg2_WD ( I, J, DEP_HG ) + CALL ADD_Hg2_SNOWPACK( I, J, DEP_HG, & State_Met, State_Chm, State_Diag ) ! Check if it is a HgP tag ELSE IF ( SpcInfo%Is_HgP ) THEN - ! Get the category # for this HgP species - Hg_Cat = SpcInfo%Hg_Cat - ! Archive wet-deposited HgP - CALL ADD_HgP_WD ( I, J, Hg_Cat, DEP_HG ) - CALL ADD_Hg2_SNOWPACK( I, J, Hg_Cat, DEP_HG, & + CALL ADD_HgP_WD ( I, J, DEP_HG ) + CALL ADD_Hg2_SNOWPACK( I, J, DEP_HG, & State_Met, State_Chm, State_Diag ) ENDIF diff --git a/Headers/species_database_mod.F90 b/Headers/species_database_mod.F90 index ed1a8604b..f7ee424d1 100644 --- a/Headers/species_database_mod.F90 +++ b/Headers/species_database_mod.F90 @@ -540,7 +540,6 @@ SUBROUTINE Init_Species_Database( Input_Opt, SpcData, SpcCount, RC ) IF ( RC /= GC_SUCCESS ) GOTO 999 IF ( v_bool ) THEN SpcCount%nHg0 = SpcCount%nHg0 + 1 - ThisSpc%Hg_Cat = SpcCount%nHg0 ThisSpc%Is_Hg0 = v_bool ENDIF @@ -549,7 +548,6 @@ SUBROUTINE Init_Species_Database( Input_Opt, SpcData, SpcCount, RC ) IF ( RC /= GC_SUCCESS ) GOTO 999 IF ( v_bool ) THEN SpcCount%nHg2 = SpcCount%nHg2 + 1 - ThisSpc%Hg_Cat = SpcCount%nHg2 ThisSpc%Is_Hg2 = v_bool ENDIF @@ -558,7 +556,6 @@ SUBROUTINE Init_Species_Database( Input_Opt, SpcData, SpcCount, RC ) IF ( RC /= GC_SUCCESS ) GOTO 999 IF ( v_bool ) THEN SpcCount%nHgP = SpcCount%nHgP + 1 - ThisSpc%Hg_Cat = SpcCount%nHgP ThisSpc%Is_HgP = v_bool ENDIF diff --git a/Headers/species_mod.F90 b/Headers/species_mod.F90 index 287bca614..740d2e87c 100644 --- a/Headers/species_mod.F90 +++ b/Headers/species_mod.F90 @@ -156,10 +156,9 @@ MODULE Species_Mod LOGICAL :: MP_SizeResNum ! T=size-resolved aerosol number ! Tagged mercury parameters - LOGICAL :: Is_Hg0 ! T=total or tagged Hg0 species - LOGICAL :: Is_Hg2 ! T=total or tagged Hg2 species - LOGICAL :: Is_HgP ! T=total or tagged HgP species - INTEGER :: Hg_Cat ! Tagged Hg category number + LOGICAL :: Is_Hg0 ! Is a Hg0 species? + LOGICAL :: Is_Hg2 ! Is a Hg2 species? + LOGICAL :: Is_HgP ! Is a HgP species? END TYPE Species ! @@ -379,7 +378,6 @@ SUBROUTINE Spc_Zero( Spc ) Spc%DryAltId = MISSING_INT Spc%DryDepId = MISSING_INT Spc%GasSpcId = MISSING_INT - Spc%Hg_Cat = MISSING_INT Spc%HygGrthId = MISSING_INT Spc%KppFixId = MISSING_INT Spc%KppSpcId = MISSING_INT diff --git a/Headers/state_chm_mod.F90 b/Headers/state_chm_mod.F90 index 934f3dec9..e17ae3dfc 100644 --- a/Headers/state_chm_mod.F90 +++ b/Headers/state_chm_mod.F90 @@ -195,22 +195,16 @@ MODULE State_Chm_Mod !----------------------------------------------------------------------- ! For the tagged Hg simulation !----------------------------------------------------------------------- - INTEGER :: N_HG_CATS ! # of Hg categories - INTEGER, POINTER :: Hg0_Id_List(: ) ! Hg0 cat <-> tracer # - INTEGER, POINTER :: Hg2_Id_List(: ) ! Hg2 cat <-> tracer # - INTEGER, POINTER :: HgP_Id_List(: ) ! HgP cat <-> tracer # - CHARACTER(LEN=4), POINTER :: Hg_Cat_Name(: ) ! Category names - - REAL(fp), POINTER :: OceanHg0(:,:,:) ! Hg(0) ocean mass [kg] - REAL(fp), POINTER :: OceanHg2(:,:,:) ! Hg(II) ocean mass [kg] - REAL(fp), POINTER :: OceanHgP(:,:,:) ! HgP ocean mass [kg] - REAL(fp), POINTER :: SnowHgOcean(:,:,:) ! Reducible Hg snowpack + REAL(fp), POINTER :: OceanHg0(:,:) ! Hg(0) ocean mass [kg] + REAL(fp), POINTER :: OceanHg2(:,:) ! Hg(II) ocean mass [kg] + REAL(fp), POINTER :: OceanHgP(:,:) ! HgP ocean mass [kg] + REAL(fp), POINTER :: SnowHgOcean(:,:) ! Reducible Hg snowpack ! on ocean [kg] - REAL(fp), POINTER :: SnowHgLand(:,:,:) ! Reducible Hg snowpack + REAL(fp), POINTER :: SnowHgLand(:,:) ! Reducible Hg snowpack ! on land [kg] - REAL(fp), POINTER :: SnowHgOceanStored(:,:,:) ! Non-reducible Hg + REAL(fp), POINTER :: SnowHgOceanStored(:,:) ! Non-reducible Hg ! snowpack on ocean - REAL(fp), POINTER :: SnowHgLandStored(:,:,:) ! Non-reducible Hg + REAL(fp), POINTER :: SnowHgLandStored(:,:) ! Non-reducible Hg ! snowpack on land !---------------------------------------------------------------------- @@ -477,11 +471,6 @@ SUBROUTINE Zero_State_Chm( State_Chm, RC ) State_Chm%IsorropBisulfate => NULL() ! Hg simulation quantities - State_Chm%N_HG_CATS = 0 - State_Chm%Hg0_Id_List => NULL() - State_Chm%Hg2_Id_List => NULL() - State_Chm%HgP_Id_List => NULL() - State_Chm%Hg_Cat_Name => NULL() State_Chm%OceanHg0 => NULL() State_Chm%OceanHg2 => NULL() State_Chm%OceanHgP => NULL() @@ -2016,7 +2005,7 @@ SUBROUTINE Init_State_Chm( Input_Opt, State_Chm, State_Grid, RC ) ENDIF !======================================================================= - ! Initialize State_Chm quantities pertinent to Hg/tagHg simulations + ! Initialize State_Chm quantities pertinent to Hg simulations !======================================================================= IF ( Input_Opt%ITS_A_MERCURY_SIM ) THEN CALL Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & @@ -2496,78 +2485,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & thisLoc = & ' -> at Init_Hg_Simulation_Fields (in module Headers/state_chm_mod.F90)' - ! Hg0, Hg2, HgP should all have the same number of categories as - ! returned from the species database. If not, there's an error. -!>> IF ( SpcCount%nHg0 == SpcCount%nHg2 .and. & -!>> SpcCount%nHg0 == SpcCount%nHgP ) THEN - State_Chm%N_Hg_CATS = 1 !SpcCount%nHg0 -!>> ELSE -!>> ErrMsg = 'Inconsistent number of Hg categories!' -!>> CALL GC_Error( ErrMsg, RC, ThisLoc ) -!>> RETURN -!>> ENDIF - - ! Index array: Hg0 species # <--> Hg0 category # - ALLOCATE( State_Chm%Hg0_Id_List( State_Chm%N_Hg_CATS ), STAT=RC ) - CALL GC_CheckVar( 'State_Chm%Hg0_Id_List', 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%Hg0_Id_List = 0 - -!>> ! Index array: Hg2 species # <--> Hg0 category # -!>> ALLOCATE( State_Chm%Hg2_Id_List( State_Chm%N_Hg_CATS ), STAT=RC ) -!>> CALL GC_CheckVar( 'State_Chm%Hg2_Id_List', 0, RC ) -!>> IF ( RC /= GC_SUCCESS ) RETURN -!>> State_Chm%Hg2_Id_List = 0 -!>> -!>> ! Index array: HgP species # <--> Hg0 category # -!>> ALLOCATE( State_Chm%HgP_Id_List( State_Chm%N_Hg_CATS ), STAT=RC ) -!>> CALL GC_CheckVar( 'State_Chm%HgP_Id_List', 0, RC ) -!>> IF ( RC /= GC_SUCCESS ) RETURN -!>> State_Chm%HgP_Id_List = 0 - - ! Hg category names - ALLOCATE( State_Chm%Hg_Cat_Name( State_Chm%N_Hg_CATS ), STAT=RC ) - CALL GC_CheckVar( 'State_Chm%Hg_Cat_Name', 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%Hg_Cat_Name = '' - - ! Loop over all species - DO N = 1, State_Chm%nSpecies - - ! Point to Species Database entry for Hg species N - ThisSpc => State_Chm%SpcData(N)%Info - - ! Populate the Hg0 index array - IF ( ThisSpc%Is_Hg0 ) THEN - State_Chm%Hg0_Id_List(ThisSpc%Hg_Cat) = ThisSpc%ModelId - ENDIF - - ! Populate the Hg2 index array - IF ( ThisSpc%Is_Hg2 ) THEN - State_Chm%Hg2_Id_List(ThisSpc%Hg_Cat) = ThisSpc%ModelId - ENDIF - - ! Populate the HgP index array - IF ( ThisSpc%Is_HgP ) THEN - State_Chm%HgP_Id_List(ThisSpc%Hg_Cat) = ThisSpc%ModelId - ENDIF - - ! Free pointer - ThisSpc => NULL() - ENDDO - - ! Loop over Hg categories (except the first) - DO C = 2, State_Chm%N_Hg_CATS - - ! Hg0 tracer number corresponding to this category - N = State_Chm%Hg0_Id_List(C) - - ! The category name (e.g. "_can") follows the "Hg0" - ThisSpc => State_Chm%SpcData(N)%Info - State_Chm%Hg_Cat_Name(C) = ThisSpc%Name(4:7) - ThisSpc => NULL() - ENDDO - !------------------------------------------------------------------------ ! Hg(0) ocean mass !------------------------------------------------------------------------ @@ -2578,7 +2495,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%OceanHg0, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2597,7 +2513,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%OceanHg2, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2616,7 +2531,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%OceanHgP, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2635,7 +2549,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%SnowHgOcean, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2654,7 +2567,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%SnowHgLand, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2673,7 +2585,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%SnowHgOceanStored, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) IF ( RC /= GC_SUCCESS ) THEN @@ -2692,7 +2603,6 @@ SUBROUTINE Init_Hg_Simulation_Fields( Input_Opt, State_Chm, State_Grid, & State_Grid = State_Grid, & chmId = chmId, & Ptr2Data = State_Chm%SnowHgLandStored, & - nSlots = State_Chm%N_Hg_CATS, & RC = RC ) !------------------------------------------------------------------------ @@ -2940,34 +2850,6 @@ SUBROUTINE Cleanup_State_Chm( State_Chm, RC ) State_Chm%BoundaryCond => NULL() ENDIF - IF ( ASSOCIATED( State_Chm%Hg_Cat_Name ) ) THEN - DEALLOCATE( State_Chm%Hg_Cat_Name, STAT=RC ) - CALL GC_CheckVar( 'State_Chm%Hg_Cat_Name', 2, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%Hg_Cat_Name => NULL() - ENDIF - - IF ( ASSOCIATED( State_Chm%Hg0_Id_List ) ) THEN - DEALLOCATE( State_Chm%Hg0_Id_List, STAT=RC ) - CALL GC_CheckVar( 'State_Chm%Hg0_Id_List', 2, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%Hg0_Id_List => NULL() - ENDIF - - IF ( ASSOCIATED( State_Chm%Hg2_Id_List ) ) THEN - DEALLOCATE( State_Chm%Hg2_Id_List, STAT=RC ) - CALL GC_CheckVar( 'State_Chm%Hg2_Id_List', 2, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%Hg2_Id_List => NULL() - ENDIF - - IF ( ASSOCIATED( State_Chm%HgP_Id_List ) ) THEN - DEALLOCATE( State_Chm%HgP_Id_List, STAT=RC ) - CALL GC_CheckVar( 'State_Chm%HgP_Id_List', 2, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Chm%HgP_Id_List => NULL() - ENDIF - IF ( ASSOCIATED( State_Chm%AeroArea ) ) THEN DEALLOCATE( State_Chm%AeroArea, STAT=RC ) CALL GC_CheckVar( 'State_Chm%AeroArea', 2, RC ) @@ -4184,43 +4066,36 @@ SUBROUTINE Get_Metadata_State_Chm( am_I_Root, metadataID, Found, & IF ( isDesc ) Desc = 'Hg(0) ocean mass' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'OCEANHG2' ) IF ( isDesc ) Desc = 'Hg(II) ocean mass' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'OCEANHGP' ) IF ( isDesc ) Desc = 'HgP ocean mass' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'SNOWHGOCEAN' ) IF ( isDesc ) Desc = 'Reducible Hg snowpack on ocean' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'SNOWHGLAND' ) IF ( isDesc ) Desc = 'Reducible Hg snowpack on land' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'SNOWHGOCEANSTORED' ) IF ( isDesc ) Desc = 'Non-reducible Hg snowpack on ocean' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE( 'SNOWHGLANDSTORED' ) IF ( isDesc ) Desc = 'Non-reducible Hg snowpack on land' IF ( isUnits ) Units = 'kg' IF ( isRank ) Rank = 2 - IF ( isSpc ) PerSpc = 'HgCat' CASE ( 'IODIDE' ) IF ( isDesc ) Desc = 'Surface iodide concentration' @@ -4952,7 +4827,7 @@ FUNCTION Test_for_Species_Dim( perSpc ) RESULT( returnCode ) !------------------------------------------------------------------------------ !BOC SELECT CASE( TRIM( perSpc ) ) - CASE( 'ADV', 'ALL', 'DRY', 'WET', 'HgCat' ) + CASE( 'ADV', 'ALL', 'DRY', 'WET' ) returnCode = 1 CASE( '' ) returnCode = 0 @@ -5002,8 +4877,6 @@ FUNCTION Get_NumSlots( perSpc, State_Chm ) RESULT( nSlots ) nSlots = State_Chm%nDryDep CASE( 'WET' ) nSlots = State_Chm%nWetDep - CASE( 'HgCat' ) - nSlots = State_Chm%N_HG_CATS CASE DEFAULT nSlots = -1 END SELECT @@ -5054,49 +4927,27 @@ SUBROUTINE Get_Diagnostic_Name( State_Chm, perSpc, N, name, & ! Objects TYPE(Species), POINTER :: ThisSpc - IF ( PerSpc == 'HgCat' ) THEN - - !--------------------------------------------------------------------- - ! Hg simulation quantities - !--------------------------------------------------------------------- - -!>> ! Append the category name to the diagnostic name -!>> diagName = TRIM( name ) // TRIM( State_Chm%Hg_Cat_Name(N) ) -!>> -!>> ! Append the category name to the description -!>> diagDesc = TRIM( desc ) // TRIM( State_Chm%Hg_Cat_Name(N) ) - - ! Append the species name to the diagnostic name with an underscore - diagName = TRIM( name )! // '_' // TRIM( State_Chm%Hg_Cat_Name(N) ) + !--------------------------------------------------------------------- + ! All other species-bound quantities + !--------------------------------------------------------------------- - ! Append the species name to the diagnostic description - diagDesc = TRIM( desc )! // ' ' // TRIM( State_Chm%Hg_Cat_Name(N) ) + ! Get the species index from the diagnostic index + ! depending on the value of PerSpc (bmy, 05 Oct 2021) + modelId = N + IF ( PerSpc == 'DRY' ) modelId = State_Chm%Map_DryDep(N) + IF ( PerSpc == 'WET' ) modelId = State_Chm%Map_WetDep(N) + + ! Point to the proper species, by modelId + ThisSpc => State_Chm%SpcData(modelId)%Info - ELSE - - !--------------------------------------------------------------------- - ! All other species-bound quantities - !--------------------------------------------------------------------- - - ! Get the species index from the diagnostic index - ! depending on the value of PerSpc (bmy, 05 Oct 2021) - modelId = N - IF ( PerSpc == 'DRY' ) modelId = State_Chm%Map_DryDep(N) - IF ( PerSpc == 'WET' ) modelId = State_Chm%Map_WetDep(N) - - ! Point to the proper species, by modelId - ThisSpc => State_Chm%SpcData(modelId)%Info + ! Append the species name to the diagnostic name with an underscore + diagName = TRIM( name ) // '_' // TRIM( ThisSpc%Name ) - ! Append the species name to the diagnostic name with an underscore - diagName = TRIM( name ) // '_' // TRIM( ThisSpc%Name ) + ! Append the species name to the diagnostic description + diagDesc = TRIM( desc ) // ' ' // TRIM( ThisSpc%Name ) - ! Append the species name to the diagnostic description - diagDesc = TRIM( desc ) // ' ' // TRIM( ThisSpc%Name ) - - ! Free pointer - ThisSpc => NULL() - - ENDIF + ! Free pointer + ThisSpc => NULL() END SUBROUTINE Get_Diagnostic_Name !EOC diff --git a/run/GCClassic/HISTORY.rc.templates/HISTORY.rc.Hg b/run/GCClassic/HISTORY.rc.templates/HISTORY.rc.Hg index cf33c46de..a4afc781b 100644 --- a/run/GCClassic/HISTORY.rc.templates/HISTORY.rc.Hg +++ b/run/GCClassic/HISTORY.rc.templates/HISTORY.rc.Hg @@ -29,6 +29,7 @@ EXPID: ./OutputDir/GEOSChem #============================================================================== COLLECTIONS: 'Restart', #'MercuryChem', + #'MercuryEmis', #'MercuryOcean', 'SpeciesConc', #'Budget', @@ -63,13 +64,13 @@ COLLECTIONS: 'Restart', 'Met_DELPDRY ', 'Met_BXHEIGHT ', 'Met_TropLev ', - #'Chem_OceanHg0 ', - #'Chem_OceanHg2 ', - #'Chem_OceanHgP ', - #'Chem_SnowHgOcean ', - #'Chem_SnowHgLand ', - #'Chem_SnowHgOceanStored ', - #'Chem_SnowHgLandStored ', + 'Chem_OceanHg0 ', + 'Chem_OceanHg2 ', + 'Chem_OceanHgP ', + 'Chem_SnowHgOcean ', + 'Chem_SnowHgLand ', + 'Chem_SnowHgOceanStored ', + 'Chem_SnowHgLandStored ', :: #============================================================================== # %%%%% THE SpeciesConc COLLECTION %%%%% @@ -187,6 +188,23 @@ COLLECTIONS: 'Restart', 'Hg2GasToSSA ', :: #============================================================================== +# %%%%% THE MercuryEmis COLLECTION %%%%% +# +# Concentration and prod/loss diagnostics for mercury species that are +# not archived by HEMCO. +# +# Available for the Hg and tagged Hg simulations +#============================================================================== + MercuryEmis.template: '%y4%m2%d2_%h2%n2z.nc4', + MercuryEmis.frequency: ${RUNDIR_HIST_TIME_AVG_FREQ} + MercuryEmis.duration: ${RUNDIR_HIST_TIME_AVG_DUR} + MercuryEmis.mode: 'time-averaged' + MercuryEmis.fields: 'EmisHg0land ', + 'EmisHg0ocean ', + 'EmisHg0snow ', + 'EmisHg0soil ', +:: +#============================================================================== # %%%%% THE MercuryOcean COLLECTION %%%%% # # Oceanic masses and fluxes of mercury species diff --git a/run/shared/species_database_hg.yml b/run/shared/species_database_hg.yml index e300efb10..a4c2aa0a9 100644 --- a/run/shared/species_database_hg.yml +++ b/run/shared/species_database_hg.yml @@ -44,59 +44,58 @@ Hg0: Is_Hg0: true MW_g: 200.59 Hg_CHEM_PROP: &HgChemProperties + DD_F0: 0.0 + DD_Hstar: 1.0e+14 Henry_CR: 8.40e+03 Henry_K0: 1.40e+06 Is_Advected: true Is_DryDep: true Is_Gas: true - Is_Photolysis: true Is_WetDep: true + Is_Hg2: true WD_RetFactor: 1.0 + MW_g: 200.59 HgBr: Fullname: HgBr Formula: HgBr Is_Advected: true Is_Gas: true Is_Photolysis: true - MW_g: 280.49 + MW_g: 200.59 HgBrNO2: + << : *HgChemProperties Fullname: syn-HgBrONO Formula: BrHgONO - Is_Advected: true - Is_DryDep: true - Is_Gas: true Is_Photolysis: true - MW_g: 326.50 HgBrHO2: << : *HgChemProperties Fullname: HgBrHO2 Formula: BrHgOOH - MW_g: 313.50 + Is_Photolysis: true HgBrBrO: << : *HgChemProperties Fullname: HgBrBrO Formula: BrHgOBr - MW_g: 376.40 + Is_Photolysis: true HgBrClO: << : *HgChemProperties Fullname: HgBrClO Formula: BrHgOCl - MW_g: 331.95 + Is_Photolysis: true HgBrO: FullName: HgBrO Formula: HgBrO Is_Gas: true - MW_g: 296.49 HgBrOH: << : *HgChemProperties Fullname: HgBrOH Formula: BrHgOH - MW_g: 297.50 + Is_Photolysis: true HgBr2: << : *HgChemProperties Fullname: HgBr2 Formula: HgBr2 - MW_g: 360.40 + Is_Photolysis: true HgCl: Fullname: HgCl Formula: HgCl @@ -104,96 +103,95 @@ HgCl: Is_Gas: true Is_Photolysis: true Is_WetDep: false - MW_g: 236.04 + MW_g: 200.59 HgClNO2: << : *HgChemProperties Fullname: syn-HgClONO Formula: ClHgONO - MW_g: 282.00 + Is_Photolysis: true HgClHO2: - WD_RetFactor: 1.0 << : *HgChemProperties Fullname: HgClHO2 Formula: ClHgOOH - MW_g: 269.00 - WD_RetFactor: 1.0 + Is_Photolysis: true HgClClO: << : *HgChemProperties Fullname: HgClClO Formula: ClHgOCl - MW_g: 287.50 + Is_Photolysis: true HgClBrO: << : *HgChemProperties Fullname: HgClBrO Formula: ClHgOBr - MW_g: 332.00 + Is_Photolysis: true HgClBr: << : *HgChemProperties Fullname: HgClBr Formula: HgBrCl - MW_g: 316.00 + Is_Photolysis: true HgClO: Fullname: HgClO Formula: ClHgO Is_Gas: true - MW_g: 252.04 + MW_g: 200.59 HgClOH: << : *HgChemProperties Fullname: HgClOH Formula: ClHgOH - MW_g: 253.00 + Is_Photolysis: true HgOH: Fullname: HgOH Formula: HgOH Is_Advected: true Is_Gas: true Is_Photolysis: true - MW_g: 201.00 + MW_g: 200.59 HgOHO: Fullname: HgOHO Formula: HgOHO Is_Gas: true - MW_g: 233.60 + MW_g: 200.59 HgOHNO2: << : *HgChemProperties Fullname: syn-HgOHONO Formula: HOHgONO - MW_g: 263.60 + Is_Photolysis: true HgOHHO2: << : *HgChemProperties Fullname: HgOHHO2 Formula: HOHgOOH - MW_g: 250.60 + Is_Photolysis: true HgOHClO: << : *HgChemProperties Fullname: HgBrClO Formula: HOHgOCl - MW_g: 269.0000 + Is_Photolysis: true HgOHBrO: << : *HgChemProperties Fullname: HgOHBrO Formula: HOHgOBr - MW_g: 313.5000 + Is_Photolysis: true HgOHOH: << : *HgChemProperties Fullname: HgOH2 Formula: HOHgOH - MW_g: 234.60 + Is_Photolysis: false HgCl2: << : *HgChemProperties Fullname: HgCl2 Formula: HgCl2 - MW_g: 271.50 + Is_Photolysis: false Hg2ClP: Fullname: Hg(II) chloride salts on sea-salt aerosols Formula: HgCln Is_Aerosol: true Is_DryDep: true Is_WetDep: true - MW_g: 201.00 + Is_HgP: true + MW_g: 200.59 WD_AerScavEff: 1.0 WD_KcScaleFac: [1.0, 0.5, 1.0] - WD_RainoutEff: [1.0, 1.0, 1.0] + WD_RainoutEff: [1.0, 0.0, 1.0] Hg2ORGP: Fullname: Hg(II) organic complex in aerosols Formula: R-Hg @@ -202,17 +200,18 @@ Hg2ORGP: Is_DryDep: true Is_Photolysis: true Is_WetDep: true - MW_g: 201.00 - WD_AerScavEff: 1.0 + Is_HgP: true + MW_g: 200.59 + WD_AerScavEff: 0.8 WD_KcScaleFac: [1.0, 0.5, 1.0] - WD_RainoutEff: [1.0, 1.0, 1.0] + WD_RainoutEff: [1.0, 0.0, 1.0] Hg2STRP: Fullname: Hg(II) in stratospheric aerosols Formula: Hg2+ Is_Advected: true Is_Aerosol: true - MW_g: 201.00 - MW_g: 35.45 + Is_HgP: true + MW_g: 200.59 HO2: Background_VV: 4.0e-15 Formula: HO2