!------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: hcox_gc_POPs_mod.F90 ! ! !DESCRIPTION: Defines the HEMCO extension for the GEOS-Chem persistent ! organic pollutants (POPs) specialty simulation. !\\ !\\ ! !INTERFACE: ! MODULE HCOX_GC_POPs_Mod ! ! !USES: ! USE HCO_Error_Mod USE HCO_Diagn_Mod USE HCO_State_Mod, ONLY : HCO_State ! Derived type for HEMCO state USE HCOX_State_Mod, ONLY : Ext_State ! Derived type for External state IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! PUBLIC :: HcoX_GC_POPs_Run PUBLIC :: HcoX_GC_POPs_Init PUBLIC :: HcoX_Gc_POPs_Final ! ! !PRIVATE MEMBER FUNCTIONS: ! PRIVATE :: VEGEMISPOP PRIVATE :: LAKEEMISPOP PRIVATE :: SOILEMISPOP PRIVATE :: IS_LAND PRIVATE :: IS_ICE ! ! !REMARKS: ! ! References: ! ============================================================================ ! (1 ) Zhang, Y., and S. Tao, Global atmospheric emission inventory of ! polycyclic aromatic hydrocarbons (PAHs) for 2004. Atm Env, 43, 812-819, ! 2009. ! (2 ) Friedman, C.L, and N.E. Selin, Long-Range Atmospheric Transport of ! Polycyclic Aromatic Hydrocarbons: A Global 3-D Model Analysis ! Including Evaluation of Arctic Sources, Environ. Sci. Technol., 46(17), ! 9501-9510, 2012. ! (3 ) Friedman, C.L., Y. Zhang, and N.E. Selin, Climate change and ! emissions impacts on atmospheric PAH transport to the Arctic, Environ. ! Sci. Technol., 48, 429-437, 2014. ! (4 ) Friedman, C.L., J.R. Pierce, and N.E. Selin, Assessing the influence of ! secondary organic versus primary carbonaceous aerosols on long-range ! atmospheric polycyclic aromatic hydrocarbon transport, Environ. Sci. ! Technol., 48(6), 3293-3302, 2014. ! ! !REVISION HISTORY: ! 20 Sep 2010 - N.E. Selin - Initial Version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !DEFINED PARAMETERS: ! REAL(hp), PARAMETER :: SMALLNUM = 1e-20_hp ! ! !PRIVATE TYPES: ! TYPE :: MyInst ! Fields required by module INTEGER :: Instance INTEGER :: ExtNr ! HEMCO Extension number INTEGER :: IDTPOPPOCPO ! Index # for POPPOC tracer INTEGER :: IDTPOPPBCPO ! Index # for POPPBC tracer INTEGER :: IDTPOPG ! Index # for POPG tracer ! Pointers to emission arrays read from disk REAL(sp), POINTER :: POP_TOT_EM(:,:) => NULL() ! [kg/m2/s] REAL(sp), POINTER :: POP_SURF(:,:) => NULL() ! [kg] REAL(sp), POINTER :: C_OC(:,:,:) => NULL() ! [kg] REAL(sp), POINTER :: C_BC(:,:,:) => NULL() ! [kg] REAL(sp), POINTER :: F_OC_SOIL(:,:) => NULL() ! [kg/m2] ! Calculated emissions of OC-phase, BC-phase, and gas-phase POPs [kg/m2/s] REAL(sp), POINTER :: EmisPOPPOCPO(:,:,:) REAL(sp), POINTER :: EmisPOPPBCPO(:,:,:) REAL(sp), POINTER :: EmisPOPG (:,:,:) ! For diagnostics REAL(sp), POINTER :: EmisPOPGfromSoil (:,:) REAL(sp), POINTER :: EmisPOPGfromLake (:,:) REAL(sp), POINTER :: EmisPOPGfromLeaf (:,:) REAL(sp), POINTER :: EmisPOPGfromSnow (:,:) REAL(sp), POINTER :: EmisPOPGfromOcean (:,:) REAL(sp), POINTER :: FluxPOPGfromSoilToAir(:,:) REAL(sp), POINTER :: FluxPOPGfromAirToSoil(:,:) REAL(sp), POINTER :: FluxPOPGfromLakeToAir(:,:) REAL(sp), POINTER :: FluxPOPGfromAirToLake(:,:) REAL(sp), POINTER :: FluxPOPGfromLeafToAir(:,:) REAL(sp), POINTER :: FluxPOPGfromAirtoLeaf(:,:) REAL(sp), POINTER :: FugacitySoilToAir (:,:) REAL(sp), POINTER :: FugacityLakeToAir (:,:) REAL(sp), POINTER :: FugacityLeafToAir (:,:) TYPE(MyInst), POINTER :: NextInst => NULL() END TYPE MyInst ! Pointer to instances TYPE(MyInst), POINTER :: AllInst => NULL() CONTAINS !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: HCOX_GC_POPs_run ! ! !DESCRIPTION: Subroutine HcoX\_Gc\_POPs\_Run computes emissions of OC-phase, ! BC-phase, and gas-phase POPs for the GEOS-Chem POPs specialty simulation. !\\ !\\ ! !INTERFACE: ! SUBROUTINE HCOX_GC_POPs_Run( ExtState, HcoState, RC ) ! ! !USES: ! USE HCO_EmisList_Mod, ONLY : HCO_GetPtr USE HCO_FluxArr_Mod, ONLY : HCO_EmisAdd USE HCO_Clock_Mod, ONLY : HcoClock_First ! ! !INPUT PARAMETERS: ! TYPE(Ext_State), POINTER :: ExtState ! Options for POPs sim TYPE(HCO_State), POINTER :: HcoState ! HEMCO state ! ! !INPUT/OUTPUT PARAMETERS: ! INTEGER, INTENT(INOUT) :: RC ! Success or failure? ! ! !REMARKS: ! This code is based on routine EMISSPOPS in prior versions of GEOS-Chem. ! ! !REVISION HISTORY: ! 20 Sep 2010 - N.E. Selin - Initial Version based on EMISSMERCURY ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !DEFINED PARAMETERS: ! ! Universal gas constant for adjusting KOA for temp: 8.3144598 [J/mol/K] REAL(hp), PARAMETER :: R = 8.3144598e+0_hp ! Density of octanol, needed for partitioning into OC: 820 [kg/m^3] REAL(hp), PARAMETER :: DENS_OCT = 82e+1_hp ! Density of BC, needed for partitioning onto BC: 1 [kg/L] or 1000 [kg/m^3] ! From Lohmann and Lammel, Environ. Sci. Technol., 2004, 38:3793-3803. REAL(hp), PARAMETER :: DENS_BC = 1e+3_hp ! ! !LOCAL VARIABLES: ! INTEGER :: I, J, L INTEGER :: PBL_MAX INTEGER :: MONTH, YEAR REAL(hp) :: F_OF_PBL, TK REAL(hp) :: T_POP REAL(hp) :: C_OC1, C_BC1 REAL(hp) :: C_OC2, C_BC2 REAL(hp) :: F_POP_OC, F_POP_BC REAL(hp) :: F_POP_G, AIR_VOL REAL(hp) :: KOA_T, KBC_T REAL(hp) :: KOC_BC_T, KBC_OC_T REAL(hp) :: VR_OC_AIR, VR_OC_BC REAL(hp) :: VR_BC_AIR, VR_BC_OC REAL(hp) :: SUM_F REAL(hp) :: OC_AIR_RATIO, OC_BC_RATIO REAL(hp) :: BC_AIR_RATIO, BC_OC_RATIO REAL(hp) :: FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND REAL(hp) :: FRAC_LEAF, FRAC_LAKE, FRAC_SOIL ! LOGICAL, SAVE :: FIRST = .TRUE. LOGICAL :: aIR LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE CHARACTER(LEN=255) :: MSG, LOC ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer ! from gas phase to OC. For now we use Delta H for phase transfer ! from the gas phase to the pure liquid state. ! For PHENANTHRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol]. ! For PYRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol]. ! For BENZO[a]PYRENE: ! this is also taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol] REAL(hp) :: DEL_H ! KOA_298 for partitioning of gas phase POP to atmospheric OC ! KOA_298 = Cpop in octanol/Cpop in atmosphere at 298 K ! For PHENANTHRENE: ! log KOA_298 = 7.64, or 4.37*10^7 [unitless] ! For PYRENE: ! log KOA_298 = 8.86, or 7.24*10^8 [unitless] ! For BENZO[a]PYRENE: ! log KOA_298 = 11.48, or 3.02*10^11 [unitless] ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825). REAL(hp) :: KOA_298 ! KBC_298 for partitioning of gas phase POP to atmospheric BC ! KBC_298 = Cpop in black carbon/Cpop in atmosphere at 298 K ! For PHENANTHRENE: ! log KBC_298 = 10.0, or 1.0*10^10 [unitless] ! For PYRENE: ! log KBC_298 = 11.0, or 1.0*10^11 [unitless] ! For BENZO[a]PYRENE: ! log KBC_298 = 13.9, or 7.94*10^13 [unitless] ! (Lohmann and Lammel, EST, 2004, 38:3793-3802) REAL(hp) :: KBC_298 ! Pointers REAL(sp), POINTER :: Arr3D(:,:,:) TYPE(MyInst), POINTER :: Inst !======================================================================= ! HCOX_GC_POPs_RUN begins here! !======================================================================= LOC = 'HCOX_GC_POPs_RUN (HCOX_GC_POPS_MOD.F90)' ! Return if extension not turned on IF ( ExtState%GC_POPs <= 0 ) RETURN ! Enter CALL HCO_ENTER( HcoState%Config%Err, LOC, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC ) RETURN ENDIF ! Get instance Inst => NULL() CALL InstGet ( ExtState%GC_POPs, Inst, RC ) IF ( RC /= HCO_SUCCESS ) THEN WRITE(MSG,*) 'Cannot find GC_POPs instance Nr. ', ExtState%GC_POPs CALL HCO_ERROR(MSG,RC) RETURN ENDIF DEL_H = ExtState%POP_DEL_H KOA_298 = ExtState%POP_KOA KBC_298 = ExtState%POP_KBC Arr3D => NULL() !======================================================================= ! Get pointers to gridded data imported through config. file !======================================================================= IF ( HcoClock_First(HcoState%Clock,.TRUE.) ) THEN CALL HCO_GetPtr( HcoState, 'TOT_POP', Inst%POP_TOT_EM, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC ) RETURN ENDIF CALL HCO_GetPtr( HcoState, 'GLOBAL_OC', Inst%C_OC, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC ) RETURN ENDIF CALL HCO_GetPtr( HcoState, 'GLOBAL_BC', Inst%C_BC, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC ) RETURN ENDIF CALL HCO_GetPtr( HcoState, 'SURF_POP', Inst%POP_SURF, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC ) RETURN ENDIF CALL HCO_GetPtr( HcoState, 'SOIL_CARBON', Inst%F_OC_SOIL, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC ) RETURN ENDIF ! Convert F_OC_SOIL from kg/m2 to fraction DO J=1, HcoState%NY DO I=1, HcoState%NX ! Assume most of carbon mass extends to 5 cm and calculate ! concentration in kg/kg ! For now, assume a mean soil bulk density of 1300 kg/m3 similar to ! McLachlan 2002 to calculate a dry weight fraction Inst%F_OC_SOIL(I,J) = Inst%F_OC_SOIL(I,J) / 30e-2_hp / 13e+2_hp ENDDO ENDDO ! FIRST = .FALSE. ENDIF ! Maximum extent of the PBL [model level] IF ( .NOT. ASSOCIATED(ExtState%PBL_MAX) ) THEN CALL HCO_ERROR ( 'PBL_MAX not defined in ExtState!', RC ) RETURN ELSE PBL_MAX = DBLE( ExtState%PBL_MAX ) ENDIF !================================================================= ! Call emissions routines for revolatilization fluxes from surfaces ! Assume all re-emited POPs are in the gas phase until partitioning ! to ambient OC and BC in the boundary layer. ! Re-emission flux/mass depends on type of POP ! First draft, CLF, 28 Aug 2012 !================================================================= CALL SOILEMISPOP( Inst%POP_SURF, Inst%F_OC_SOIL, HcoState, ExtState, Inst ) CALL LAKEEMISPOP( Inst%POP_SURF, HcoState, ExtState, Inst ) CALL VEGEMISPOP ( Inst%POP_SURF, HcoState, ExtState, Inst ) Inst%EmisPOPGfromSnow = 0e+0_hp Inst%EmisPOPGfromOcean = 0e+0_hp ! Loop over grid boxes DO J = 1, HcoState%Ny DO I = 1, HcoState%Nx F_OF_PBL = 0e+0_hp T_POP = 0e+0_hp ! Here, save the total from the emissions array ! into the T_POP variable [kg/m2/s] T_POP = Inst%POP_TOT_EM(I,J) ! Now add revolatilization (secondary) emissions to primary [kg/m2/s] T_POP = T_POP + Inst%EmisPOPGfromSoil(I,J) + & Inst%EmisPOPGfromLake(I,J) + & Inst%EmisPOPGfromLeaf(I,J) !+ & !Inst%EmisPOPGfromSnow(I,J) + & !Inst%EmisPOPGfromOcean(I,J) !==================================================================== ! Apportion total POPs emitted to gas phase, OC-bound, and BC-bound ! emissions (clf, 2/1/2011) ! Then partition POP throughout PBL; store into STT [kg] ! Now make sure STT does not underflow (cdh, bmy, 4/6/06; eck 9/20/10) !==================================================================== ! Loop up to max PBL level DO L = 1, PBL_MAX ! Get temp [K] TK = ExtState%TK%Arr%Val(I,J,L) ! Define temperature-dependent partition coefficients: ! KOA_T, the octanol-air partition coeff at temp T [unitless] KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp))) ! Define KBC_T, the BC-air partition coeff at temp T [unitless] ! TURN OFF TEMPERATURE DEPENDENCY FOR SENSITIVITY ANALYSIS KBC_T = KBC_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp))) ! Define KOC_BC_T, theoretical OC-BC part coeff at temp T [unitless] KOC_BC_T = KOA_T / KBC_T ! Define KBC_OC_T, theoretical BC_OC part coeff at temp T [unitless] KBC_OC_T = 1d0 / KOC_BC_T ! Get monthly mean OC and BC concentrations [kg/box] C_OC1 = Inst%C_OC(I,J,L) C_BC1 = Inst%C_BC(I,J,L) ! Make sure OC is not negative C_OC1 = MAX( C_OC1, 0e+0_hp ) ! Convert C_OC and C_BC units to volume per box ! [m^3 OC or BC/box] C_OC2 = C_OC1 / DENS_OCT C_BC2 = C_BC1 / DENS_BC ! Get air volume (m^3) AIR_VOL = ExtState%AIRVOL%Arr%Val(I,J,L) ! Define volume ratios: ! VR_OC_AIR = volume ratio of OC to air [unitless] VR_OC_AIR = C_OC2 / AIR_VOL ! VR_OC_BC = volume ratio of OC to BC [unitless] VR_OC_BC = C_OC2 / C_BC2 ! VR_BC_AIR = volume ratio of BC to air [unitless] VR_BC_AIR = VR_OC_AIR / VR_OC_BC ! VR_BC_OC = volume ratio of BC to OC [unitless] !VR_BC_OC(I,J,L) = 1d0 / VR_OC_BC(I,J,L) VR_BC_OC = 1d0 / VR_OC_BC ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase, ! BC-phase, and gas phase with new time step (should only change if ! temp changes or OC/BC concentrations change) OC_AIR_RATIO = 1e+0_hp / (KOA_T * VR_OC_AIR) OC_BC_RATIO = 1e+0_hp / (KOC_BC_T * VR_OC_BC) BC_AIR_RATIO = 1e+0_hp / (KBC_T * VR_BC_AIR) BC_OC_RATIO = 1e+0_hp / (KBC_OC_T * VR_BC_OC) ! If there are zeros in OC or BC concentrations, make sure they ! don't cause problems with phase fractions IF ( C_OC1 > SMALLNUM .and. C_BC1 > SMALLNUM ) THEN F_POP_OC = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO + OC_BC_RATIO) F_POP_BC = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO + BC_OC_RATIO) ELSE IF ( C_OC1 > SMALLNUM .and. C_BC1 .le. SMALLNUM ) THEN F_POP_OC = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO) F_POP_BC = SMALLNUM ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 > SMALLNUM ) THEN F_POP_OC = SMALLNUM F_POP_BC = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO) ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 .le. SMALLNUM) THEN F_POP_OC = SMALLNUM F_POP_BC = SMALLNUM ENDIF ! Gas-phase: F_POP_G = 1e+0_hp - F_POP_OC - F_POP_BC ! Check that sum of fractions equals 1 SUM_F = F_POP_OC + F_POP_BC + F_POP_G ! Fraction of PBL that box (I,J,L) makes up [unitless] F_OF_PBL = ExtState%FRAC_OF_PBL%Arr%Val(I,J,L) ! Calculate rates of POP emissions in each phase [kg/m2/s] ! OC-phase: Inst%EmisPOPPOCPO(I,J,L) = F_POP_OC * F_OF_PBL * T_POP ! BC-phase Inst%EmisPOPPBCPO(I,J,L) = F_POP_BC * F_OF_PBL * T_POP ! Gas-phase Inst%EmisPOPG(I,J,L) = F_POP_G * F_OF_PBL * T_POP ENDDO ENDDO ENDDO !======================================================================= ! Add POPs emissions to HEMCO data structure & diagnostics !======================================================================= !---------------------- ! OC-PHASE EMISSIONS !---------------------- IF ( Inst%IDTPOPPOCPO > 0 ) THEN ! Add flux to emissions array Arr3D => Inst%EmisPOPPOCPO(:,:,:) CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPOCPO, & RC, ExtNr=Inst%ExtNr ) Arr3D => NULL() IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPOCPO', RC ) RETURN ENDIF ENDIF !---------------------- ! BC-PHASE EMISSIONS !---------------------- IF ( Inst%IDTPOPPBCPO > 0 ) THEN ! Add flux to emissions array Arr3D => Inst%EmisPOPPBCPO(:,:,:) CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPBCPO, & RC, ExtNr=Inst%ExtNr ) Arr3D => NULL() IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPBCPO', RC ) RETURN ENDIF ENDIF !---------------------- ! GASEOUS EMISSIONS !---------------------- IF ( Inst%IDTPOPG > 0 ) THEN ! Add flux to emissions array Arr3D => Inst%EmisPOPG(:,:,:) CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPG, & RC, ExtNr=Inst%ExtNr ) Arr3D => NULL() IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPG', RC ) RETURN ENDIF ENDIF !======================================================================= ! Cleanup & quit !======================================================================= ! Nullify pointers Inst => NULL() ! Return w/ success CALL HCO_LEAVE( HcoState%Config%Err,RC ) END SUBROUTINE HCOX_GC_POPs_Run !EOC !------------------------------------------------------------------------------ ! GEOS-Chem Global Chemical Transport Model ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: soilemispop ! ! !DESCRIPTION: Subroutine SOILEMISPOP is the subroutine for secondary ! POP emissions from soils. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SOILEMISPOP( POP_SURF, F_OC_SOIL, HcoState, ExtState, Inst ) ! ! !INPUT PARAMETERS: ! REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg] REAL(sp), DIMENSION(:,:), INTENT(IN) :: F_OC_SOIL ! Frac C in soil [g/g] TYPE(HCO_STATE), POINTER :: HcoState ! Hemco state TYPE(Ext_State), POINTER :: ExtState ! Module options TYPE(MyInst), POINTER :: Inst ! Instance ! ! !REVISION HISTORY: ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! LOGICAL :: IS_SNOW_OR_ICE LOGICAL :: IS_LAND_OR_ICE INTEGER :: I, J, L REAL(hp) :: POPG REAL(hp) :: TK_SURF REAL(hp) :: KSA_T, FLUX, F_OC REAL(hp) :: SOIL_CONC, DTSRCE, KOA_T REAL(hp) :: DIFF, FSOIL, FAIR, DS REAL(hp) :: DSA, DAD, DWD, PL REAL(hp) :: ZSOIL, ZAIR, TK REAL(hp) :: DTCHEM, NEWSOIL, E_KDEG REAL(hp) :: FRAC_LAKE, FRAC_SOIL REAL(hp) :: FUG_R REAL(hp) :: AREA_M2 ! ! !DEFINED PARAMETERS: ! ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer ! from gas phase to OC. For now we use Delta H for phase transfer ! from the gas phase to the pure liquid state. ! For PHENANTHRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol]. ! For PYRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol]. ! For BENZO[a]PYRENE: ! this is also taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol] REAL(hp) :: DEL_H ! R = universal gas constant for adjusting KOA for temp: ! 8.3144598 [J/mol/K OR m3*Pa/K/mol] REAL(hp), PARAMETER :: R = 8.3144598d0 ! Molecular weight ! For phe, 0.17823 kg/mol ! For pyr, 0.20225 kg/mol ! For BaP, 0,25231 kg/mol REAL(hp) :: MW ! Molecular weight of air REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol ! For PHENANTHRENE: ! log KOA_298 = 7.64, or 4.37*10^7 [unitless] ! For PYRENE: ! log KOA_298 = 8.86, or 7.24*10^8 [unitless] ! For BENZO[a]PYRENE: ! log KOA_298 = 11.48, or 3.02*10^11 [unitless] ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825). REAL(hp) :: KOA_298 ! Set transfer velocity and diffusion coefficient values REAL(hp), PARAMETER :: KSA = 1d0 ![m/h] REAL(hp), PARAMETER :: BA = 0.04d0 ![m2/h] REAL(hp), PARAMETER :: BW = 4d-6 ![m2/h] ! Set soil degradation rate REAL(hp), PARAMETER :: DEGR = 3.5d-5 ![/h] !================================================================= ! SOILEMISPOP begins here! !================================================================= ! Copy values from ExtState DEL_H = ExtState%POP_DEL_H KOA_298 = ExtState%POP_KOA MW = ExtState%POP_XMW ! Emission timestep [s] DTSRCE = HcoState%TS_EMIS ! Chemistry timestep [h] DTCHEM = HcoState%TS_CHEM / 60e+0_hp DO J=1, HcoState%NY DO I=1, HcoState%NX ! Set logicals ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE) ! IS_LAND will return non-ocean boxes but may still contain lakes ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE) IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR. & ( IS_ICE (I,J,ExtState) ) ) IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR. & ( IS_LAND(I,J,ExtState) .AND. & ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) ) ! Do soils routine only if we are on land that is not covered with ! snow or ice IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN ! Get fraction of grid box covered by lake surface area FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J) ! Get fraction of land remaining ! Assume the remaining land is soil and get OC content. ! If remaining land is not soil (e.g., desert), there ! should be a characteristically low OC content ! that will have little capacity to store POPs ! ONLY SUBTRACT FRAC LAKE NOW FRAC_SOIL = MAX(1e+0_hp - FRAC_LAKE, 0e+0_hp) ! Get surface skin temp [K] TK_SURF = ExtState%TSKIN%Arr%Val(I,J) ! Get air temp [K] TK = ExtState%TK%Arr%Val(I,J,1) ! Get gas phase air POP concentration at surface in mol/m3 ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15) POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM ) ! (kg trc/kg dry air) / (0.178 kg trc/mol) * (kg dry air/m3) POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3 ! Grid box surface area [m2] AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J) ! Get soil concentration in top 5 cm of soil ! (following Howsam et al 2000) ! From Howsam et al, soil burdens are equal to ! 2.6 years deposition for PHE, ! 10 years for PYR, and 9.4 years for BaP ! Convert to mol/m3 ! 2.6 yrs * kg deposited to soil in 1 yr * / 0.178 kg/mol ! / area grid box (m2) / 0.05 m SOIL_CONC = 10e+0_hp * POP_SURF(I,J) / MW / & AREA_M2 / 5e-2_hp ! mol/m3 ! Get rid of mass due to degradation ! Use rate constant for BaP from Mackay and Paterson 1991: ! 3.5*10^-5 /h ! Calculate exponential factor E_KDEG = EXP (-DEGR * DTCHEM) ! Adjust conc NEWSOIL = SOIL_CONC * E_KDEG ! Get foc from GTMM saved files F_OC = F_OC_SOIL(I,J) ! Define temperature-dependent KOA: KOA_T = KOA_298 * EXP((-DEL_H/R) * ( ( 1e+0_hp / TK_SURF ) - & ( 1e+0_hp / 298e+0_hp ) )) ! Dimensionless coefficient (mol/m3 soil / mol/m3 air) ! KSA = 1.5 (fTOC)*Koa KSA_T = 1.5 * F_OC * KOA_T KSA_T = MAX( KSA_T, SMALLNUM ) ! Calculate fugacities from concentrations by dividing by "Z" ! values, or the fugacity capacity in mol/m3*Pa following ! Mackay and Paterson, 1991 ! fsoil = Csoil * R * T / KSA [Pa] ! where Csoil is in mol/m3, R is in Pa * m3 / mol K ! T is in K and KSA is dimensionless FSOIL = NEWSOIL * R * TK_SURF / KSA_T ! fair = Cair * R * T [Pa] ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K FAIR = POPG * R * TK ! Calculate the fugacity gradient [Pa] ! If the gradient is negative, fair is larger and the POP will ! diffuse from air to soil ! If the gradient is positive, fsoil is larger and POP will ! diffuse from soil to air DIFF = FSOIL - FAIR FUG_R = FSOIL/FAIR ! Calculate "Z" values from fugacities. ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f ZAIR = POPG / FAIR ! (mol/m3) / (Pa) ZSOIL = NEWSOIL / FSOIL ! (mol/m3) / (Pa) ! Calculate the "D" value, or the transfer coefficient that ! describes the movement of POP between phases (Mackay and ! Paterson, 1991). [mol/h*Pa] ! The D value for soil-air diffusion is given by ! Ds = 1 / (1/Dsa + 1/(Dad + Dwd)) ! Dsa is the air-side boundary layer diffusion parameter [mol/h*Pa] ! Dad is the diffusion parameter between soil particles and ! "soil air" [mol/h*Pa] ! Dwd is the diffusion parameter between soil particles and ! porewater [mol/h*Pa] ! Dsa is in series with soil-air and soil-water diffusion, ! which are in parallel ! Need to define each D value ! DSA = kSA * Zair ! where kSA is a mass transfer coefficient [m/h], ! Zair is the air fugacity capacity [mol/m3*Pa] ! *********** ! DAD = BA * Zair ! where BA is the molecular diffusivity in air [m2/h] ! *********** ! DWD = BW * Zsoil ! where BW is the molecular diffusivity in water [m2/h] ! **** PL = the soil diffusion pathlength, set to half the ! soil depth (0.025 m) DSA = KSA * ZAIR ! (m/h) * (mol/m3*Pa) = mol/m2*h*Pa DAD = BA* ZAIR ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa DWD = BW * ZSOIL ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa PL = 0.025e+0_hp DS = 1e+0_hp / ( 1e+0_hp/DSA + PL/(DAD+DWD) ) ! mol/(m2*h*Pa) [* m3*Pa/K/mol = m/h/K ! Calculate Flux in mol/m2/h FLUX = DS * DIFF ! Change to units of ng/m2/d for storage FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp ! Kludge soil emissions from poles for now ! Bug somewhere that allows GCAP versions to think some high polar ! boxes during some months are land rather than ice - results in ! extremely high fluxes IF ( HcoState%Grid%YMID%Val(I,J) > 60 .OR. & HcoState%Grid%YMID%Val(I,J) < -60 ) THEN FLUX = 0e+0_hp DIFF = 0e+0_hp ENDIF ! Convert to an emission rate in kg/m2/s for returning to ! HcoX_GC_POPs_Run Inst%EmisPOPGfromSoil(I,J) = MAX( FLUX / 24e+0_hp / 3600e+0_hp / & 1e+12_hp, 0e+0_hp ) ! Multiply the emissions by the fraction of land that is soil Inst%EmisPOPGfromSoil(I,J) = Inst%EmisPOPGfromSoil(I,J) * FRAC_SOIL ! If the flux is positive, then the direction will be from the ! soil to the air. ! If the flux is zero or negative, store it in a separate array. IF ( FLUX > 0e+0_hp ) THEN ! Store positive flux Inst%FluxPOPGfromSoilToAir(I,J) = FLUX ! Make sure negative flux diagnostic has nothing added to it Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp ! Store the soil/air fugacity ratio Inst%FugacitySoilToAir(I,J) = FUG_R ELSE IF ( FLUX <= 0e+0_hp ) THEN ! Store the negative flux Inst%FluxPOPGfromAirToSoil(I,J) = FLUX ! Add nothing to positive flux Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp ! Continue to store the fugacity ratio Inst%FugacitySoilToAir(I,J) = FUG_R ENDIF ELSE ! We are not on land or the land is covered with ice or snow FLUX = 0e+0_hp Inst%EmisPOPGfromSoil(I,J) = 0e+0_hp Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp Inst%FugacitySoilToAir(I,J) = 0e+0_hp ENDIF ENDDO ENDDO END SUBROUTINE SOILEMISPOP !EOC !------------------------------------------------------------------------------ ! GEOS-Chem Global Chemical Transport Model ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: lakeemispop ! ! !DESCRIPTION: Subroutine LAKEEMISPOP is the subroutine for secondary ! POP emissions from lakes. !\\ !\\ ! !INTERFACE: ! SUBROUTINE LAKEEMISPOP( POP_SURF, HcoState, ExtState, Inst ) ! ! !INPUT PARAMETERS: ! REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg] TYPE(HCO_STATE), POINTER :: HcoState ! Hemco state TYPE(Ext_State), POINTER :: ExtState ! Module options TYPE(MyInst), POINTER :: Inst ! Instance ! ! !REVISION HISTORY: ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! INTEGER :: I, J, L REAL(hp) :: TK_SURF REAL(hp) :: KAW_T, FLUX, KOL_T REAL(hp) :: DTSRCE REAL(hp) :: KA_H2O, KA_POP, KW_CO2, KW_POP REAL(hp) :: DA_H2O, DA_POP, DW_CO2, DW_POP REAL(hp) :: SCH_CO2, SCH_POP REAL(hp) :: TK, PRESS REAL(hp) :: C_DISS REAL(hp) :: POPG, U10M, ALPHA REAL(hp) :: FRAC_LAKE, VISC_H2O REAL(hp) :: SFCWINDSQR REAL(hp) :: AREA_M2 LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE ! ! !DEFINED PARAMETERS: ! ! Delta H for POP: ! For PHENANTHRENE: ! this is the Delta H for phase transfer ! from air to water (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47 [kJ/mol]. ! For PYRENE: 43000 [kJ/mol] REAL(hp) :: DEL_HW ! R = universal gas constant for adjusting KOA for temp: ! 8.3144598d-3 [kJ/mol/K] REAL(hp), PARAMETER :: R = 8.3144598d-3 ! Molecular weight ! For phe, 0.17823 kg/mol REAL(hp) :: MWPOP ! Molecular weight of air REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol ! Molecular weight of water REAL(hp), PARAMETER :: MWH2O = 18.1d0 ! g/mol ! Molar volumes calculated following Abraham and McGowan 1987 as ! summarized by Schwarzenbach et al. 2003. ! Each element is assigned a characteristic atomic volume, and an atomic ! volume of 6.56 cm3/mol is subtracted for each bond, no matter whether ! single, double, or triple ! C = 16.35, H = 8.71, O = 12.43, N = 14.39, P = 24.87, F = 10.48 ! Br = 26.21, I = 34.53, S = 22.91, Si = 26.83 ! Molar volume of water ! 2*(8.71) + 12.43 - 2*(6.56) = 16.73 REAL(hp), PARAMETER :: V_H2O = 16.73d0 ! cm3/mol ! Molar volume of CO2 ! 16.35 + 2*(12.43) - 2*(6.56) = 28.1 REAL(hp), PARAMETER :: V_CO2 = 28.1d0 ! cm3/mol ! Molar volume of POP ! For PHE (C16H10): ! 16*(16.35) + 10*(8.71) - 29*(6.56) = REAL(hp), PARAMETER :: V_POP = 538.94d0 ! cm3/mol ! Molar volume of air - average of gases in air REAL(hp), PARAMETER :: V_AIR = 20.1d0 ! cm3/mol ! For PHENANTHRENE: ! log KAW_298 = -2.76, or 1.74*10-3 [unitless] ! For PYRENE: ! log KAW_298 = -3.27, or 5.37*10-4 [unitless] REAL(hp) :: KAW_298 ! Set the kinematic viscosity of freshwater at 20C !REAL(hp), PARAMETER :: VISC_H2O ! = 1d0 ![cm2/s] ! Set aqueous degradation rate !REAL(hp), PARAMETER :: DEGR != 3.5d-5 ![/h] !================================================================= ! LAKEEMISPOP begins here! !================================================================= ! Copy values from ExtState DEL_HW = ExtState%POP_DEL_HW MWPOP = ExtState%POP_XMW KAW_298 = ExtState%POP_HSTAR ! Emission timestep [s] DTSRCE = HcoState%TS_EMIS DO J=1, HcoState%NY DO I=1, HcoState%NX ! Set logicals ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE) ! IS_LAND will return non-ocean boxes but may still contain lakes ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE) IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR. & ( IS_ICE (I,J,ExtState) ) ) IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR. & ( IS_LAND(I,J,ExtState) .AND. & ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) ) ! Do soils routine only if we are on land that is not covered with ! snow or ice IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN ! Get fraction of grid box covered by lake surface area FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J) IF ( FRAC_LAKE > 0e+0_hp ) THEN ! Get surface skin temp [K] TK_SURF = ExtState%TSKIN%Arr%Val(I,J) ! Get air temp [K] TK = ExtState%TK%Arr%Val(I,J,1) ! Get surface pressure at end of dynamic time step [hPa] PRESS = ExtState%PSC2_WET%Arr%Val(I,J) ! Convert to units of atm PRESS = PRESS / 1013.25e+0_hp ! Get gas phase air POP concentration at surface in mol/m3 ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15) POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM ) ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3) POPG = POPG / MWPOP * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3 ! Grid box surface area [m2] AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J) ! Get the dissolved POP concentration at lake surface in mol/m3 ! Distribute the total deposited mass to a volume that best ! matches observed dissolved concentrations ! Future versions should consider aqueous particle concentrations ! and sinking rates and photolytic/microbial degradation ! Start with 1 m - scale by 100 C_DISS = POP_SURF(I,J) / MWPOP / AREA_M2 / 100e+0_hp ! Wind speed at 10m altitude [m/s] SFCWINDSQR = ExtState%U10M%Arr%Val(I,J)**2 & + ExtState%V10M%Arr%Val(I,J)**2 U10M = SQRT( SFCWINDSQR ) ! Need to calculate water-side and air-side mass transfer ! coefficients ! Start with air-side ! Relate air-side MTC of POP to that of H2O ! First, calculate air-side MTC for water following ! Schwarzenbach, Gschwend, Imboden 2003 KA_H2O = 0.2e+0_hp * U10M + 0.3e+0_hp ! cm/s ! Relate air-side MTC for water to that of POP via diffusivities ! following Schwarzenbach et al 2003 ! Calculate temperature-dependent diffusivities in air first ! folling Fuller et al. 1966 (summarized by Schwarzenbach ! et al. 2003) DA_POP = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + & (1e+0_hp / (MWPOP*1e+3_hp) ) )**0.5e+0_hp / & ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) + & V_POP**(1e+0_hp/3e+0_hp))**2e+0_hp ) ! cm2/s DA_H2O = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + & (1e+0_hp / MWH2O ) )**0.5e+0_hp / & ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) + & V_H2O**(1e+0_hp/3e+0_hp))**2e+0_hp ) ! cm2/s ! Relate POP and H2O air-side MTCs KA_POP = KA_H2O * ( DA_POP / DA_H2O )**( 0.67e+0_hp ) ! cm/s ! Now calculate water-side MTCs ! Start with calculating the water side MTC of CO2 ! This depends on wind speed (Schwarzenbach et al 2003) ! Three different scenarios for the water surface under ! different wind speeds are considered: ! Smooth Surface Regime (SSR): u10 <= 4.2 m/s ! Rough Surface Regime (RSR): 4.2 m/s < u10 <= 13 m/s ! Breaking Wave Regime (BWR): u10 > 13 m/s ! Alpha, the exponent in the relationship for the water side MTC, ! also depends on the wind speed. Set this as well IF ( U10M <= 4.2e+0_hp ) THEN KW_CO2 = 0.65e-3_hp ! cm/s ALPHA = 0.67e+0_hp ELSE IF ( U10M > 4.2e+0_hp .AND. U10M <= 13e+0_hp ) THEN KW_CO2 = ( 0.79e+0_hp * U10M - 2.68e+0_hp ) * 1e-3_hp ALPHA = 0.5e+0_hp ELSE IF (U10M > 13e+0_hp) THEN KW_CO2 = ( 1.64e+0_hp * U10M - 13.69e+0_hp ) * 1e-3_hp ALPHA = 0.50e+0_hp ENDIF ! Get the temperature-dependent kinematic viscosity of water IF (TK_SURF <= 273.15 ) THEN VISC_H2O = 1.787e-2_hp ! [cm2/s] ELSE IF (TK_SURF > 273.15 .AND. TK_SURF <= 278.15 ) THEN VISC_H2O = 1.518e-2_hp ELSE IF (TK_SURF > 258.15 .AND. TK_SURF <= 283.15 ) THEN VISC_H2O = 1.307e-2_hp ELSE IF (TK_SURF > 283.15 .AND. TK_SURF <= 287.15 ) THEN VISC_H2O = 1.139e-2_hp ELSE IF (TK_SURF > 287.15 .AND. TK_SURF <= 293.15 ) THEN VISC_H2O = 1.002e-2_hp ELSE IF (TK_SURF > 293.15 .AND. TK_SURF <= 298.15 ) THEN VISC_H2O = 0.89e-2_hp ELSE IF (TK_SURF > 298.15 ) THEN VISC_H2O = 0.797e-2_hp ENDIF ! Calculate the diffusivites of CO2 and POP in water DW_CO2 = ( 13.26 * 1e-5_hp ) / & (( VISC_H2O*1e+2_hp )**1.14e+0_hp * & (V_CO2)**0.589e+0_hp ) ! [cm2/s] DW_POP = ( 13.26 * 1e-5_hp ) / & (( VISC_H2O*1e+2_hp )**1.14e+0_hp * & (V_POP)**0.589e+0_hp ) ! [cm2/s] ! Calculate the Schmidt numbers for CO2 and POP SCH_CO2 = VISC_H2O / DW_CO2 ! [unitless] SCH_POP = VISC_H2O / DW_POP ! [unitless] ! Calculate the water-side MTC for POP KW_POP = KW_CO2 * ( SCH_POP / SCH_CO2 ) ** (-ALPHA) ! [cm/s] ! Calculate the temperature-dependent dimensionless Henry's Law ! constant KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - & (1e+0_hp/298e+0_hp))) ! [unitless] ! Now calculate the overall air-water MTC KOL_T = 1e+0_hp / ( 1e+0_hp/KW_POP + & 1e+0_hp / (KA_POP*KAW_T) ) ! [cm/s] ! Calculate Flux in ng/m2/day ! FLUX = KOL_T * 3600e+0_hp * 24e+0_hp * ( C_DISS - & POPG/KAW_T ) * MWPOP * 1e+12_hp / 100e+0_hp ! Convert to an emission rate in kg/m2/s for returning to ! HcoX_GC_POPs_Run ! Only return it if it's positive Inst%EmisPOPGfromLake(I,J) = MAX(FLUX / 24e+0_hp / 3600e+0_hp / & 1e+12_hp, 0e+0_hp ) ! Multiply the emissions by the fraction of land that is water Inst%EmisPOPGfromLake(I,J) = Inst%EmisPOPGfromLake(I,J) * & FRAC_LAKE ! If the flux is positive, then the direction will be from the ! soil to the air. ! If the flux is zero or negative, store it in a separate array. IF ( FLUX > 0e+0_hp ) THEN ! Store positive flux Inst%FluxPOPGfromLakeToAir(I,J) = + FLUX ! Make sure negative flux diagnostic has nothing added to it Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp ! Store the soil/air fugacity ratio Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T) ELSE IF ( FLUX <= 0e+0_hp ) THEN ! Store the negative flux Inst%FluxPOPGfromAirToLake(I,J) = FLUX ! Add nothing to positive flux Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp ! Continue to store the fugacity ratio Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T) ENDIF ENDIF ELSE ! We are not on land or the land is covered with ice or snow ! or we are land but there is no water FLUX = 0e+0_hp Inst%EmisPOPGfromLake = 0e+0_hp Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp Inst%FugacityLakeToAir(I,J) = 0e+0_hp ENDIF ENDDO ENDDO END SUBROUTINE LAKEEMISPOP !EOC !------------------------------------------------------------------------------ ! GEOS-Chem Global Chemical Transport Model ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: vegemispop ! ! !DESCRIPTION: Subroutine VEGEMISPOP is the subroutine for secondary ! POP emissions from leaves. !\\ !\\ ! !INTERFACE: ! SUBROUTINE VEGEMISPOP( POP_SURF, HcoState, ExtState, Inst ) ! ! !INPUT PARAMETERS: ! REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg] TYPE(HCO_State), POINTER :: HcoState ! Hemco state TYPE(Ext_State), POINTER :: ExtState ! Module options TYPE(MyInst), POINTER :: Inst ! Instance ! ! !REVISION HISTORY: ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! INTEGER :: I, J, L REAL(hp) :: POPG, POPG_GL REAL(hp) :: FRAC_SNOWFREE_LAND, TK_SURF REAL(hp) :: KLA_T, FLUX, K_MT REAL(hp) :: LEAF_CONC, DTSRCE, KOA_T REAL(hp) :: DIFF, FLEAF, FAIR, DS REAL(hp) :: DAB_F, DC, DAL, PC, UC REAL(hp) :: ZLEAF, ZAIR, TK REAL(hp) :: NEWLEAF REAL(hp) :: LAI, KAW_T, KOW_T REAL(hp) :: FUG_R, NEW_LEAF, DLA REAL(hp) :: AREA_M2 LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE LOGICAL :: IS_SNOWFREE_LAND ! ! !DEFINED PARAMETERS: ! ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer ! from gas phase to OC. For now we use Delta H for phase transfer ! from the gas phase to the pure liquid state. ! For PHENANTHRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol]. ! For PYRENE: ! this is taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol]. ! For BENZO[a]PYRENE: ! this is also taken as the negative of the Delta H for phase transfer ! from the pure liquid state to the gas phase (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol] REAL(hp) :: DEL_H ! Delta H for POP [kJ/mol]. ! For PHENANTHRENE: ! this is the Delta H for phase transfer ! from air to water (Schwarzenbach, ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47000 [J/mol]. ! For PYRENE: 43000 [J/mol] REAL(hp) :: DEL_HW ! R = universal gas constant for adjusting KOA for temp: ! 8.314 [J/mol/K OR m3*Pa/K/mol] REAL(hp), PARAMETER :: R = 8.3144598e+0_hp ! Molecular weight ! For phe, 0.17823 kg/mol REAL(hp) :: MW ! Molecular weight of air REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol ! For PHENANTHRENE: ! log KOA_298 = 7.64, or 4.37*10^7 [unitless] ! For PYRENE: ! log KOA_298 = 8.86, or 7.24*10^8 [unitless] ! For BENZO[a]PYRENE: ! log KOA_298 = 11.48, or 3.02*10^11 [unitless] ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825). REAL(hp) :: KOA_298 ! For PHENANTHRENE: ! log KAW_298 = -2.76, or 1.74*10-3 [unitless] ! For PYRENE: ! log KAW_298 = -3.27, or 5.37*10-4 [unitless] REAL(hp) :: KAW_298 ! Set volume fractions of octanol and water in surface and reservoir ! leaf compartments [unitless] REAL(hp), PARAMETER :: OCT_SURF = 0.8e+0_hp REAL(hp), PARAMETER :: OCT_RES = 0.02e+0_hp REAL(hp), PARAMETER :: H2O_RES = 0.7e+0_hp ! Set thickness of different leaf compartments. Volumes calculated by ! multiplying thicknesses by leaf area index REAL(hp), PARAMETER :: SURF_THICK = 2e-6_hp ! m REAL(hp), PARAMETER :: RES_THICK = 250e-6_hp ! m ! Set transfer velocity and diffusion coefficient values REAL(hp), PARAMETER :: UAB_F = 9e+0_hp ![m/h] ! Set soil degradation rate REAL(hp), PARAMETER :: DEGR = 3.5e-5_hp ![/h] !================================================================= ! VEGEMISPOP begins here! !================================================================= ! Copy values from ExtState DEL_H = ExtState%POP_DEL_H KOA_298 = ExtState%POP_KOA DEL_HW = ExtState%POP_DEL_Hw KAW_298 = ExtState%POP_HSTAR MW = ExtState%POP_XMW ! Emission timestep [s] DTSRCE = HcoState%TS_EMIS DO J=1, HcoState%NY DO I=1, HcoState%NX ! Set logicals ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE) ! IS_LAND will return non-ocean boxes but may still contain lakes ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE) IS_LAND_OR_ICE = ( (IS_LAND(I,J,ExtState)) .OR. & (IS_ICE (I,J,ExtState)) ) IS_SNOW_OR_ICE = ( (IS_ICE (I,J,ExtState)) .OR. & (IS_LAND(I,J,ExtState) .AND. & ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) ) ! Do soils routine only if we are on land that is not covered with ! snow or ice IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN ! Get fraction of grid box covered by leaf surface area ! Do not consider different vegetation types for now LAI = ExtState%LAI%Arr%Val(I,J) IF ( LAI > 0e+0_hp ) THEN ! Get surface skin temp [K] TK_SURF = ExtState%TSKIN%Arr%Val(I,J) ! Get air temp [K] TK = ExtState%TK%Arr%Val(I,J,1) ! Get gas phase air POP concentration at surface in mol/m3 ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15) POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM ) ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3) POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3 ! Grid box surface area [m2] AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J) ! Only consider partitioning into leaf surface (not reservoir) ! for now following Mackay et al 2006 Environ Sci & Pollut Res ! Include reservoir when land-atm models become dynamic ! Assume that all leaf surfaces contain an average lipid content ! of 80% (Mackay et al 2006) ! Get leaf concentration ! Convert to mol/m3 ! kg deposited to leaf in 1 yr * / 0.178 kg/mol ! / area grid box (m2) / surface thickness m LEAF_CONC = POP_SURF(I,J) / MW / AREA_M2 / SURF_THICK ! mol/m3 ! Check concentration in leaves by assuming a density similar ! to water (1 g/cm3) ! No degradation/metabolism for now. Just scale leaf ! concentrations to match flux observations NEWLEAF = LEAF_CONC/1e+4_hp !SCALING FACTOR ! Define temperature-dependent KOA: KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK_SURF) - & (1e+0_hp/298e+0_hp))) ! Calculate the temperature-dependent dimensionless Henry's Law ! constant KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - & (1e+0_hp/298e+0_hp))) ! [unitless] ! Estimate the temperature-dependent dimensionless octanol-water ! constant KOW_T = KOA_T * KAW_T ! Define dimensionless leaf surface-air partition coefficient ! (mol/m3 leaf / mol/m3 air) ! KLA = foct_surf * Koa KLA_T = OCT_SURF * KOA_T ! KLA_T = MAX( KLA_T, SMALLNUM ) ! Calculate fugacities from concentrations by dividing by "Z" ! values, or the fugacity capacity in mol/m3*Pa following Mackay ! and Paterson, 1991 ! fleaf = Cleaf * R * T / KSA [Pa] ! where Cleaf is in mol/m3, R is in Pa * m3 / mol K ! T is in K and KLA is dimensionless FLEAF = NEWLEAF * R * TK_SURF / KLA_T ! fair = Cair * R * T [Pa] ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K FAIR = POPG * R * TK ! Calculate the fugacity gradient [Pa] ! If the gradient is negative, fair is larger and the POP will ! diffuse from air to soil ! If the gradient is positive, fsoil is larger and POP will ! diffuse from soil to air DIFF = FLEAF - FAIR FUG_R = FLEAF/FAIR ! Calculate "Z" values from fugacities. ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f ZAIR = POPG / FAIR ! (mol/m3) / (Pa) ZLEAF = NEWLEAF / FLEAF ! (mol/m3) / (Pa) ! Calculate the "D" value, or the transfer coefficient that ! describes the movement of POP between phases (Mackay and ! Paterson, 1991, Cousins and Mackay 2000, internal report). ! [mol/h*Pa] ! The D value for leaf surface-air gas diffusion is given by ! Dla = 1 / (1/Dc + 1/(Dab-f)) [mol/(Pa*h)] ! where Dab-f is the boundary layer diffusion [mol/h*Pa] ! given by Dab-f = As * L * Uab-f * Za ! where As is the area of the land surface [m2], L is the leaf ! area index [m2/m2], ! Uab-f is a mass transfer coefficient for surface-air boundary ! layer diffusion [m/h], ! and Za is the fugacity capacity of the air [mol/(m3*Pa)] ! Dc is the cuticle diffusion, given by ! Dc = As * L * Uc * Zf ! where As and L are as above, Uc is the cuticle mass transfer ! coefficient [m/h], ! and Zf is the fugacity capacity of the leaf surface ! (mol/(m3*Pa)) ! Uc is given by ! Uc = 3600 * Pc * 1/Kaw ! where Pc is the cuticle permeance (m/s) and Kaw is the ! dimensionless air-water partition coefficient. ! Pc is given by ! Log Pc = ((0.704 * log Kow - 11.2) + ! (-3.47 - 2.79 * logMW + 0.970 log Kow)) / 2 ! (an average of two equations) ! Need to define each D value ! DAB_F: ! m/h * mol/(m3*Pa) = (mol/h*Pa*m2) DAB_F = UAB_F * ZAIR ! mol/(h*Pa*m2) ! Calculate PC and then Uc in order to calculate Dc ! PC, UC are calculated according to Cousins and Mackay, ! Chemosphere, 2001, Table 2 PC = 10** (( 0.704e+0_hp * LOG (KOW_T) - 11.2e+0_hp ) + & ( -3.47e+0_hp -2.79e+0_hp* LOG(MW*1000d0) + 0.97e+0_hp & * LOG(KOW_T))/2e+0_hp) ![m/s] UC = 3600e+0_hp * PC * 1e+0_hp/KAW_T ! [m/h] DC = UC * ZLEAF ! mol/(h*Pa*m2) ! Now calculate overall transfer ! mol/(h*Pa*m2) DLA = 1e+0_hp / (1e+0_hp/DC + 1e+0_hp/DAB_F) ! Calculate Flux in mol/h/m2 FLUX = DLA * DIFF ! Change to units of ng/m2/d for storage FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp ! Convert to an emission rate in kg/m2/s for returning to ! HcoX_GC_POPs_Run ! Only want to add rates that are positive Inst%EmisPOPGfromLeaf(I,J) = MAX(FLUX * LAI / 24e+0_hp / & 3600e+0_hp / 1e+12_hp, 0e+0_hp) ! If the flux is positive, then the direction will be from the ! soil to the air. ! If the flux is zero or negative, store it in a separate array. IF ( FLUX > 0e+0_hp ) THEN ! Store positive flux Inst%FluxPOPGfromLeafToAir(I,J) = FLUX ! Make sure negative flux diagnostic has nothing added to it Inst%FluxPOPGfromAirtoLeaf = 0e+0_hp ! Store the soil/air fugacity ratio Inst%FugacityLeafToAir(I,J) = FUG_R ELSE IF ( FLUX <= 0e+0_hp ) THEN ! Store the negative flux Inst%FluxPOPGfromAirtoLeaf(I,J) = FLUX ! Add nothing to positive flux Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp ! Continue to store the fugacity ratio Inst%FugacityLeafToAir(I,J) = FUG_R ENDIF ELSE ! We are not on land or the land is covered with ice or snow FLUX = 0e+0_hp Inst%EmisPOPGfromLeaf(I,J) = 0e+0_hp Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp Inst%FluxPOPGfromAirtoLeaf(I,J) = 0e+0_hp Inst%FugacityLeafToAir(I,J) = 0e+0_hp ENDIF ENDIF ENDDO ENDDO END SUBROUTINE VEGEMISPOP !EOC !------------------------------------------------------------------------------ ! GEOS-Chem Global Chemical Transport Model ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: is_land ! ! !DESCRIPTION: Function IS\_LAND returns TRUE if surface grid box (I,J) is ! a land box. !\\ !\\ ! !INTERFACE: ! FUNCTION IS_LAND( I, J, ExtState ) RESULT ( LAND ) ! ! !USES: ! USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE ! ! !INPUT PARAMETERS: ! INTEGER, INTENT(IN) :: I ! Longitude index of grid box INTEGER, INTENT(IN) :: J ! Latitude index of grid box TYPE(Ext_State), POINTER :: ExtState ! Module options ! ! !RETURN VALUE: ! LOGICAL :: LAND ! =T if it is a land box ! ! !REVISION HISTORY: ! 26 Jun 2000 - R. Yantosca - Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC LAND = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J), & ExtState%FRLANDIC%Arr%Val(I,J), & ExtState%FROCEAN%Arr%Val(I,J), & ExtState%FRSEAICE%Arr%Val(I,J), & ExtState%FRLAKE%Arr%Val(I,J) ) == 1 END FUNCTION IS_LAND !EOC !------------------------------------------------------------------------------ ! GEOS-Chem Global Chemical Transport Model ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: is_ice ! ! !DESCRIPTION: Function IS\_ICE returns TRUE if surface grid box (I,J) ! contains either land-ice or sea-ice. !\\ !\\ ! !INTERFACE: ! FUNCTION IS_ICE( I, J, ExtState ) RESULT ( ICE ) ! ! !USES: ! USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE ! ! !INPUT PARAMETERS: ! INTEGER, INTENT(IN) :: I ! Longitude index of grid box INTEGER, INTENT(IN) :: J ! Latitude index of grid box TYPE(Ext_State), POINTER :: ExtState ! Module options ! ! !RETURN VALUE: ! LOGICAL :: ICE ! =T if this is an ice box ! ! !REVISION HISTORY: ! 09 Aug 2005 - R. Yantosca - Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ICE = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J), & ExtState%FRLANDIC%Arr%Val(I,J), & ExtState%FROCEAN%Arr%Val(I,J), & ExtState%FRSEAICE%Arr%Val(I,J), & ExtState%FRLAKE%Arr%Val(I,J) ) == 2 END FUNCTION IS_ICE !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: HCOX_GC_POPs_Init ! ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Init initializes the HEMCO ! GC\_POPs extension. !\\ !\\ ! !INTERFACE: ! SUBROUTINE HCOX_GC_POPs_Init( HcoState, ExtName, ExtState, RC ) ! ! !USES: ! USE HCO_ExtList_Mod, ONLY : GetExtNr USE HCO_STATE_MOD, ONLY : HCO_GetExtHcoID ! ! !INPUT PARAMETERS: ! CHARACTER(LEN=*), INTENT(IN ) :: ExtName ! Extension name TYPE(Ext_State), POINTER :: ExtState ! Module options TYPE(HCO_State), POINTER :: HcoState ! Hemco state ! ! !INPUT/OUTPUT PARAMETERS: ! INTEGER, INTENT(INOUT) :: RC ! !REVISION HISTORY: ! 19 Aug 2014 - M. Sulprizio- Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! ! Scalars INTEGER :: I, N, nSpc, ExtNr CHARACTER(LEN=255) :: MSG, LOC CHARACTER(LEN=255) :: DiagnName CHARACTER(LEN=255) :: OutUnit ! Arrays INTEGER, ALLOCATABLE :: HcoIDs(:) CHARACTER(LEN=31), ALLOCATABLE :: SpcNames(:) ! Pointers TYPE(MyInst), POINTER :: Inst REAL(sp), POINTER :: Target2D(:,:) !======================================================================= ! HCOX_GC_POPs_INIT begins here! !======================================================================= LOC = 'HCOX_GC_POPs_INIT (HCOX_GC_POPS_MOD.F90)' ! Get the extension number ExtNr = GetExtNr( HcoState%Config%ExtList, TRIM(ExtName) ) IF ( ExtNr <= 0 ) RETURN ! Enter HEMCO CALL HCO_ENTER( HcoState%Config%Err, LOC, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC ) RETURN ENDIF ! Create Instance Inst => NULL() CALL InstCreate ( ExtNr, ExtState%GC_POPs, Inst, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR ( 'Cannot create GC_POPs instance', RC ) RETURN ENDIF ! Also fill ExtNrSS - this is the same as the parent ExtNr Inst%ExtNr = ExtNr ! Set species IDs CALL HCO_GetExtHcoID( HcoState, Inst%ExtNr, HcoIDs, SpcNames, nSpc, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC ) RETURN ENDIF ! Verbose mode IF ( HcoState%amIRoot ) THEN MSG = 'Use GC_POPs emissions module (extension module)' CALL HCO_MSG(HcoState%Config%Err,MSG ) MSG = 'Use the following species (Name: HcoID):' CALL HCO_MSG(HcoState%Config%Err,MSG) DO N = 1, nSpc WRITE(MSG,*) TRIM(SpcNames(N)), ':', HcoIDs(N) CALL HCO_MSG(HcoState%Config%Err,MSG) ENDDO ENDIF ! Set up tracer indices DO N = 1, nSpc SELECT CASE( TRIM( SpcNames(N) ) ) CASE( 'POPG', 'POPG_BaP', 'POPG_PHE', 'POPG_PYR' ) Inst%IDTPOPG = HcoIDs(N) CASE( 'POPPOCPO', 'POPPOCPO_BaP', 'POPPOCPO_PHE', 'POPPOCPO_PYR' ) Inst%IDTPOPPOCPO = HcoIDs(N) CASE( 'POPPBCPO', 'POPPBCPO_BaP', 'POPPBCPO_PHE', 'POPPBCPO_PYR' ) Inst%IDTPOPPBCPO = HcoIDs(N) CASE DEFAULT ! Do nothing END SELECT ENDDO ! ERROR: POPG tracer is not found! IF ( Inst%IDTPOPG <= 0 ) THEN RC = HCO_FAIL CALL HCO_ERROR( 'Cannot find POPG tracer in list of species!', RC ) RETURN ENDIF ! ERROR! POPPOCPO tracer is not found IF ( Inst%IDTPOPPOCPO <= 0 ) THEN RC = HCO_FAIL CALL HCO_ERROR( 'Cannot find POPPOCPO tracer in list of species!', RC ) RETURN ENDIF ! ERROR! POPPBCPO tracer is not found IF ( Inst%IDTPOPPBCPO <= 0 ) THEN RC = HCO_FAIL CALL HCO_ERROR( 'Cannot find POPPBCPO tracer in list of species!', RC ) RETURN ENDIF !======================================================================= ! Activate this module and the fields of ExtState that it uses !======================================================================= ! Activate met fields required by this extension ExtState%POPG%DoUse = .TRUE. ExtState%AIRVOL%DoUse = .TRUE. ExtState%AIRDEN%DoUse = .TRUE. ExtState%FRAC_OF_PBL%DoUse = .TRUE. ExtState%LAI%DoUse = .TRUE. ExtState%PSC2_WET%DoUse = .TRUE. ExtState%SNOWHGT%DoUse = .TRUE. ExtState%TK%DoUse = .TRUE. ExtState%TSKIN%DoUse = .TRUE. ExtState%U10M%DoUse = .TRUE. ExtState%V10M%DoUse = .TRUE. ExtState%FRLAND%DoUse = .TRUE. ExtState%FRLANDIC%DoUse = .TRUE. ExtState%FROCEAN%DoUse = .TRUE. ExtState%FRSEAICE%DoUse = .TRUE. ExtState%FRLAKE%DoUse = .TRUE. !======================================================================= ! Initialize data arrays !======================================================================= ALLOCATE( Inst%EmisPOPG( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPG', RC ) RETURN ENDIF Inst%EmisPOPG = 0.0e0_hp ALLOCATE( Inst%EmisPOPPOCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPPOCPO', RC ) RETURN ENDIF Inst%EmisPOPPOCPO = 0.0e0_hp ALLOCATE( Inst%EmisPOPPBCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPPBCPO', RC ) RETURN ENDIF Inst%EmisPOPPBCPO = 0.0e0_hp ALLOCATE( Inst%EmisPOPGfromSoil( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSoil', RC ) RETURN ENDIF Inst%EmisPOPGfromSoil = 0.0e0_hp ALLOCATE( Inst%EmisPOPGfromLake( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLake', RC ) RETURN ENDIF Inst%EmisPOPGfromLake = 0.0e0_hp ALLOCATE( Inst%EmisPOPGfromLeaf( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLeaf', RC ) RETURN ENDIF Inst%EmisPOPGfromLeaf = 0.0e0_hp ALLOCATE( Inst%EmisPOPGfromSnow( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSnow', RC ) RETURN ENDIF Inst%EmisPOPGfromSnow = 0.0e0_hp ALLOCATE( Inst%EmisPOPGfromOcean( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromOcean', RC ) RETURN ENDIF Inst%EmisPOPGfromOcean = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromSoilToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromSoilToAir', RC ) RETURN ENDIF Inst%FluxPOPGfromSoilToAir = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromAirToSoil( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToSoil', RC ) RETURN ENDIF Inst%FluxPOPGfromAirToSoil = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLakeToAir', RC ) RETURN ENDIF Inst%FluxPOPGfromLakeToAir = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromAirToLake( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLake', RC ) RETURN ENDIF Inst%FluxPOPGfromAirToLake = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLeafToAir', RC ) RETURN ENDIF Inst%FluxPOPGfromLeafToAir = 0.0e0_hp ALLOCATE( Inst%FluxPOPGfromAirToLeaf( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLeaf', RC ) RETURN ENDIF Inst%FluxPOPGfromAirToLeaf = 0.0e0_hp ALLOCATE( Inst%FugacitySoilToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FugacitySoilToAir', RC ) RETURN ENDIF Inst%FugacitySoilToAir = 0.0e0_hp ALLOCATE( Inst%FugacityLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FugacityLakeToAir', RC ) RETURN ENDIF Inst%FugacityLakeToAir = 0.0e0_hp ALLOCATE( Inst%FugacityLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC ) IF ( RC /= 0 ) THEN CALL HCO_ERROR ( 'Cannot allocate FugacityLeafToAir', RC ) RETURN ENDIF Inst%FugacityLeafToAir = 0.0e0_hp !======================================================================= ! Create manual diagnostics !======================================================================= DO I = 1,12 ! Define diagnostic names. These have to match the names ! in module HEMCO/Extensions/hcox_gc_POPs_mod.F90. IF ( I == 1 ) THEN DiagnName = 'EmisPOPGfromSoil' OutUnit = 'kg/m2/s' Target2D => Inst%EmisPOPGfromSoil ELSEIF ( I == 2 ) THEN DiagnName = 'EmisPOPGfromLake' OutUnit = 'kg/m2/s' Target2D => Inst%EmisPOPGfromLake ELSEIF ( I == 3 ) THEN DiagnName = 'EmisPOPGfromLeaf' OutUnit = 'kg/m2/s' Target2D => Inst%EmisPOPGfromLeaf ELSEIF ( I == 4 ) THEN DiagnName = 'FluxPOPGfromSoilToAir' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromSoilToAir ELSEIF ( I == 5 ) THEN DiagnName = 'FluxPOPGfromAirToSoil' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromAirToSoil ELSEIF ( I == 6 ) THEN DiagnName = 'FluxPOPGfromLakeToAir' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromLakeToAir ELSEIF ( I == 7 ) THEN DiagnName = 'FluxPOPGfromAirToLake' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromAirToLake ELSEIF ( I == 8 ) THEN DiagnName = 'FluxPOPGfromLeafToAir' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromLeafToAir ELSEIF ( I == 9 ) THEN DiagnName = 'FluxPOPGfromAirToLeaf' OutUnit = 'ng/m2/day' Target2D => Inst%FluxPOPGfromAirToLeaf ELSEIF ( I == 10 ) THEN DiagnName = 'FugacitySoilToAir' OutUnit = '1' Target2D => Inst%FugacitySoilToAir ELSEIF ( I == 11 ) THEN DiagnName = 'FugacityLakeToAir' OutUnit = '1' Target2D => Inst%FugacityLakeToAir ELSEIF ( I == 12 ) THEN DiagnName = 'FugacityLeafToAir' OutUnit = '1' Target2D => Inst%FugacityLeafToAir ENDIF ! Create manual diagnostics CALL Diagn_Create( HcoState = HcoState, & cName = TRIM( DiagnName ), & ExtNr = ExtNr, & Cat = -1, & Hier = -1, & HcoID = -1, & SpaceDim = 2, & OutUnit = OutUnit, & AutoFill = 0, & Trgt2D = Target2D, & RC = RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC ) RETURN ENDIF Target2D => NULL() ENDDO !======================================================================= ! Leave w/ success !======================================================================= IF ( ALLOCATED( HcoIDs ) ) DEALLOCATE( HcoIDs ) IF ( ALLOCATED( SpcNames ) ) DEALLOCATE( SpcNames ) ! Nullify pointers Inst => NULL() CALL HCO_LEAVE( HcoState%Config%Err,RC ) END SUBROUTINE HCOX_GC_POPs_Init !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: HCOX_GC_POPs_Final ! ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Final finalizes the HEMCO ! extension for the GEOS-Chem POPs specialty simulation. All module ! arrays will be deallocated. !\\ !\\ ! !INTERFACE: ! SUBROUTINE HCOX_GC_POPs_Final( ExtState ) ! ! !INPUT PARAMETERS: ! TYPE(Ext_State), POINTER :: ExtState ! Module options ! ! !REVISION HISTORY: ! 19 Aug 2014 - M. Sulprizio- Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC !======================================================================= ! HCOX_GC_POPs_FINAL begins here! !======================================================================= CALL InstRemove( ExtState%GC_POPs ) END SUBROUTINE HCOX_GC_POPs_Final !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: InstGet ! ! !DESCRIPTION: Subroutine InstGet returns a poiner to the desired instance. !\\ !\\ ! !INTERFACE: ! SUBROUTINE InstGet ( Instance, Inst, RC, PrevInst ) ! ! !INPUT PARAMETERS: ! INTEGER :: Instance TYPE(MyInst), POINTER :: Inst INTEGER :: RC TYPE(MyInst), POINTER, OPTIONAL :: PrevInst ! ! !REVISION HISTORY: ! 18 Feb 2016 - C. Keller - Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC TYPE(MyInst), POINTER :: PrvInst !================================================================= ! InstGet begins here! !================================================================= ! Get instance. Also archive previous instance. PrvInst => NULL() Inst => AllInst DO WHILE ( ASSOCIATED(Inst) ) IF ( Inst%Instance == Instance ) EXIT PrvInst => Inst Inst => Inst%NextInst END DO IF ( .NOT. ASSOCIATED( Inst ) ) THEN RC = HCO_FAIL RETURN ENDIF ! Pass output arguments IF ( PRESENT(PrevInst) ) PrevInst => PrvInst ! Cleanup & Return PrvInst => NULL() RC = HCO_SUCCESS END SUBROUTINE InstGet !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: InstCreate ! ! !DESCRIPTION: Subroutine InstCreate creates a new instance. !\\ !\\ ! !INTERFACE: ! SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC ) ! ! !INPUT PARAMETERS: ! INTEGER, INTENT(IN) :: ExtNr ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT( OUT) :: Instance TYPE(MyInst), POINTER :: Inst ! ! !INPUT/OUTPUT PARAMETERS: ! INTEGER, INTENT(INOUT) :: RC ! ! !REVISION HISTORY: ! 18 Feb 2016 - C. Keller - Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC TYPE(MyInst), POINTER :: TmpInst INTEGER :: nnInst !================================================================= ! InstCreate begins here! !================================================================= ! ---------------------------------------------------------------- ! Generic instance initialization ! ---------------------------------------------------------------- ! Initialize Inst => NULL() ! Get number of already existing instances TmpInst => AllInst nnInst = 0 DO WHILE ( ASSOCIATED(TmpInst) ) nnInst = nnInst + 1 TmpInst => TmpInst%NextInst END DO ! Create new instance ALLOCATE(Inst) Inst%Instance = nnInst + 1 Inst%ExtNr = ExtNr ! Attach to instance list Inst%NextInst => AllInst AllInst => Inst ! Update output instance Instance = Inst%Instance ! ---------------------------------------------------------------- ! Type specific initialization statements follow below ! ---------------------------------------------------------------- ! Return w/ success RC = HCO_SUCCESS END SUBROUTINE InstCreate !EOC !------------------------------------------------------------------------------ ! Harmonized Emissions Component (HEMCO) ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: InstRemove ! ! !DESCRIPTION: Subroutine InstRemove creates a new instance. !\\ !\\ ! !INTERFACE: ! SUBROUTINE InstRemove ( Instance ) ! ! !INPUT PARAMETERS: ! INTEGER :: Instance ! ! !REVISION HISTORY: ! 18 Feb 2016 - C. Keller - Initial version ! See https://github.com/geoschem/hemco for complete history !EOP !------------------------------------------------------------------------------ !BOC INTEGER :: RC TYPE(MyInst), POINTER :: PrevInst TYPE(MyInst), POINTER :: Inst !================================================================= ! InstRemove begins here! !================================================================= ! Init PrevInst => NULL() Inst => NULL() ! Get instance. Also archive previous instance. CALL InstGet ( Instance, Inst, RC, PrevInst=PrevInst ) ! Instance-specific deallocation IF ( ASSOCIATED(Inst) ) THEN !--------------------------------------------------------------------- ! Deallocate fields of Inst before popping off from the list ! in order to avoid memory leaks (Bob Yantosca (17 Aug 2022) !--------------------------------------------------------------------- IF ( ASSOCIATED( Inst%EmisPOPPOCPO ) ) THEN DEALLOCATE( Inst%EmisPOPPOCPO ) ENDIF Inst%EmisPOPPOCPO => NULL() IF ( ASSOCIATED( Inst%EmisPOPPBCPO ) ) THEN DEALLOCATE( Inst%EmisPOPPBCPO ) ENDIF Inst%EmisPOPPBCPO => NULL() IF ( ASSOCIATED( Inst%EmisPOPG ) ) THEN DEALLOCATE( Inst%EmisPOPG ) ENDIF Inst%EmisPOPG => NULL() IF ( ASSOCIATED( Inst%EmisPOPGfromSoil ) ) THEN DEALLOCATE( Inst%EmisPOPGfromSoil ) ENDIF Inst%EmisPOPGfromSoil => NULL() IF ( ASSOCIATED( Inst%EmisPOPGfromLake ) ) THEN DEALLOCATE( Inst%EmisPOPGfromLake ) ENDIF Inst%EmisPOPGfromLake => NULL() IF ( ASSOCIATED( Inst%EmisPOPGfromLeaf ) ) THEN DEALLOCATE( Inst%EmisPOPGfromLeaf ) ENDIF Inst%EmisPOPGfromLeaf => NULL() IF ( ASSOCIATED( Inst%EmisPOPGfromSnow ) ) THEN DEALLOCATE( Inst%EmisPOPGfromSnow ) ENDIF Inst%EmisPOPGfromSnow => NULL() IF ( ASSOCIATED( Inst%EmisPOPGfromOcean ) ) THEN DEALLOCATE( Inst%EmisPOPGfromOcean ) ENDIF Inst%EmisPOPGfromOcean => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromSoilToAir ) ) THEN DEALLOCATE( Inst%FluxPOPGfromSoilToAir ) ENDIF Inst%FluxPOPGfromSoilToAir => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromAirToSoil ) ) THEN DEALLOCATE( Inst%FluxPOPGfromAirToSoil ) ENDIF Inst%FluxPOPGfromAirToSoil => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromLakeToAir ) ) THEN DEALLOCATE( Inst%FluxPOPGfromLakeToAir ) ENDIF Inst%FluxPOPGfromLakeToAir => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromAirToLake ) ) THEN DEALLOCATE( Inst%FluxPOPGfromAirToLake ) ENDIF Inst%FluxPOPGfromAirToLake => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromLeafToAir ) ) THEN DEALLOCATE( Inst%FluxPOPGfromLeafToAir ) ENDIF Inst%FluxPOPGfromLeafToAir => NULL() IF ( ASSOCIATED( Inst%FluxPOPGfromAirtoLeaf ) ) THEN DEALLOCATE( Inst%FluxPOPGfromAirtoLeaf ) ENDIF Inst%FluxPOPGfromAirtoLeaf => NULL() IF ( ASSOCIATED( Inst%FugacitySoilToAir ) ) THEN DEALLOCATE( Inst%FugacitySoilToAir ) ENDIF Inst%FugacitySoilToAir => NULL() IF ( ASSOCIATED( Inst%FugacityLakeToAir ) ) THEN DEALLOCATE( Inst%FugacityLakeToAir ) ENDIF Inst%FugacityLakeToAir => NULL() IF ( ASSOCIATED( Inst%FugacityLeafToAir ) ) THEN DEALLOCATE( Inst%FugacityLeafToAir ) ENDIF Inst%FugacityLeafToAir => NULL() !--------------------------------------------------------------------- ! Pop off instance from list !--------------------------------------------------------------------- IF ( ASSOCIATED(PrevInst) ) THEN PrevInst%NextInst => Inst%NextInst ELSE AllInst => Inst%NextInst ENDIF DEALLOCATE(Inst) ENDIF ! Free pointers before exiting PrevInst => NULL() Inst => NULL() END SUBROUTINE InstRemove !EOC END MODULE HCOX_GC_POPs_Mod