LCOV - code coverage report
Current view: | top level - hemco/HEMCO/src/Extensions - hcox_gc_POPs_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 610 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 13 | 0.0 %

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hcox_gc_POPs_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Defines the HEMCO extension for the GEOS-Chem persistent
           9             : !   organic pollutants (POPs) specialty simulation.
          10             : !\\
          11             : !\\
          12             : ! !INTERFACE:
          13             : !
          14             : MODULE HCOX_GC_POPs_Mod
          15             : !
          16             : ! !USES:
          17             : !
          18             :   USE HCO_Error_Mod
          19             :   USE HCO_Diagn_Mod
          20             :   USE HCO_State_Mod,  ONLY : HCO_State   ! Derived type for HEMCO state
          21             :   USE HCOX_State_Mod, ONLY : Ext_State   ! Derived type for External state
          22             :
          23             :   IMPLICIT NONE
          24             :   PRIVATE
          25             : !
          26             : ! !PUBLIC MEMBER FUNCTIONS:
          27             : !
          28             :   PUBLIC  :: HcoX_GC_POPs_Run
          29             :   PUBLIC  :: HcoX_GC_POPs_Init
          30             :   PUBLIC  :: HcoX_Gc_POPs_Final
          31             : !
          32             : ! !PRIVATE MEMBER FUNCTIONS:
          33             : !
          34             :   PRIVATE :: VEGEMISPOP
          35             :   PRIVATE :: LAKEEMISPOP
          36             :   PRIVATE :: SOILEMISPOP
          37             :   PRIVATE :: IS_LAND
          38             :   PRIVATE :: IS_ICE
          39             : !
          40             : ! !REMARKS:
          41             : !
          42             : !  References:
          43             : !  ============================================================================
          44             : !  (1 ) Zhang, Y., and S. Tao, Global atmospheric emission inventory of
          45             : !        polycyclic aromatic hydrocarbons (PAHs) for 2004. Atm Env, 43, 812-819,
          46             : !        2009.
          47             : !  (2 ) Friedman, C.L, and N.E. Selin, Long-Range Atmospheric Transport of
          48             : !        Polycyclic Aromatic Hydrocarbons: A Global 3-D Model Analysis
          49             : !        Including Evaluation of Arctic Sources, Environ. Sci. Technol., 46(17),
          50             : !        9501-9510, 2012.
          51             : !  (3 ) Friedman, C.L., Y. Zhang, and N.E. Selin, Climate change and
          52             : !        emissions impacts on atmospheric PAH transport to the Arctic, Environ.
          53             : !        Sci. Technol., 48, 429-437, 2014.
          54             : !  (4 ) Friedman, C.L., J.R. Pierce, and N.E. Selin, Assessing the influence of
          55             : !        secondary organic versus primary carbonaceous aerosols on long-range
          56             : !        atmospheric polycyclic aromatic hydrocarbon transport, Environ. Sci.
          57             : !        Technol., 48(6), 3293-3302, 2014.
          58             : !
          59             : ! !REVISION HISTORY:
          60             : !  20 Sep 2010 - N.E. Selin    - Initial Version
          61             : !  See https://github.com/geoschem/hemco for complete history
          62             : !EOP
          63             : !------------------------------------------------------------------------------
          64             : !BOC
          65             : !
          66             : ! !DEFINED PARAMETERS:
          67             : !
          68             :   REAL(hp), PARAMETER           :: SMALLNUM    = 1e-20_hp
          69             : !
          70             : ! !PRIVATE TYPES:
          71             : !
          72             :   TYPE :: MyInst
          73             :    ! Fields required by module
          74             :    INTEGER               :: Instance
          75             :    INTEGER               :: ExtNr            ! HEMCO Extension number
          76             :    INTEGER               :: IDTPOPPOCPO      ! Index # for POPPOC tracer
          77             :    INTEGER               :: IDTPOPPBCPO      ! Index # for POPPBC tracer
          78             :    INTEGER               :: IDTPOPG          ! Index # for POPG   tracer
          79             :
          80             :    ! Pointers to emission arrays read from disk
          81             :    REAL(sp), POINTER     :: POP_TOT_EM(:,:) => NULL() ! [kg/m2/s]
          82             :    REAL(sp), POINTER     :: POP_SURF(:,:)   => NULL() ! [kg]
          83             :    REAL(sp), POINTER     :: C_OC(:,:,:)     => NULL() ! [kg]
          84             :    REAL(sp), POINTER     :: C_BC(:,:,:)     => NULL() ! [kg]
          85             :    REAL(sp), POINTER     :: F_OC_SOIL(:,:)  => NULL() ! [kg/m2]
          86             :
          87             :    ! Calculated emissions of OC-phase, BC-phase, and gas-phase POPs [kg/m2/s]
          88             :    REAL(sp), POINTER     :: EmisPOPPOCPO(:,:,:)
          89             :    REAL(sp), POINTER     :: EmisPOPPBCPO(:,:,:)
          90             :    REAL(sp), POINTER     :: EmisPOPG    (:,:,:)
          91             :
          92             :    ! For diagnostics
          93             :    REAL(sp), POINTER     :: EmisPOPGfromSoil     (:,:)
          94             :    REAL(sp), POINTER     :: EmisPOPGfromLake     (:,:)
          95             :    REAL(sp), POINTER     :: EmisPOPGfromLeaf     (:,:)
          96             :    REAL(sp), POINTER     :: EmisPOPGfromSnow     (:,:)
          97             :    REAL(sp), POINTER     :: EmisPOPGfromOcean    (:,:)
          98             :    REAL(sp), POINTER     :: FluxPOPGfromSoilToAir(:,:)
          99             :    REAL(sp), POINTER     :: FluxPOPGfromAirToSoil(:,:)
         100             :    REAL(sp), POINTER     :: FluxPOPGfromLakeToAir(:,:)
         101             :    REAL(sp), POINTER     :: FluxPOPGfromAirToLake(:,:)
         102             :    REAL(sp), POINTER     :: FluxPOPGfromLeafToAir(:,:)
         103             :    REAL(sp), POINTER     :: FluxPOPGfromAirtoLeaf(:,:)
         104             :    REAL(sp), POINTER     :: FugacitySoilToAir    (:,:)
         105             :    REAL(sp), POINTER     :: FugacityLakeToAir    (:,:)
         106             :    REAL(sp), POINTER     :: FugacityLeafToAir    (:,:)
         107             :
         108             :    TYPE(MyInst), POINTER :: NextInst => NULL()
         109             :   END TYPE MyInst
         110             :
         111             :   ! Pointer to instances
         112             :   TYPE(MyInst), POINTER  :: AllInst => NULL()
         113             :
         114             : CONTAINS
         115             : !EOC
         116             : !------------------------------------------------------------------------------
         117             : !                   Harmonized Emissions Component (HEMCO)                    !
         118             : !------------------------------------------------------------------------------
         119             : !BOP
         120             : !
         121             : ! !IROUTINE: HCOX_GC_POPs_run
         122             : !
         123             : ! !DESCRIPTION: Subroutine HcoX\_Gc\_POPs\_Run computes emissions of OC-phase,
         124             : !  BC-phase, and gas-phase POPs for the GEOS-Chem POPs specialty simulation.
         125             : !\\
         126             : !\\
         127             : ! !INTERFACE:
         128             : !
         129           0 :   SUBROUTINE HCOX_GC_POPs_Run( ExtState, HcoState, RC )
         130             : !
         131             : ! !USES:
         132             : !
         133             :     USE HCO_EmisList_Mod,  ONLY : HCO_GetPtr
         134             :     USE HCO_FluxArr_Mod,   ONLY : HCO_EmisAdd
         135             :     USE HCO_Clock_Mod,     ONLY : HcoClock_First
         136             : !
         137             : ! !INPUT PARAMETERS:
         138             : !
         139             :     TYPE(Ext_State),  POINTER       :: ExtState    ! Options for POPs sim
         140             :     TYPE(HCO_State),  POINTER       :: HcoState    ! HEMCO state
         141             : !
         142             : ! !INPUT/OUTPUT PARAMETERS:
         143             : !
         144             :     INTEGER,          INTENT(INOUT) :: RC          ! Success or failure?
         145             : !
         146             : ! !REMARKS:
         147             : !  This code is based on routine EMISSPOPS in prior versions of GEOS-Chem.
         148             : !
         149             : ! !REVISION HISTORY:
         150             : !  20 Sep 2010 - N.E. Selin  - Initial Version based on EMISSMERCURY
         151             : !  See https://github.com/geoschem/hemco for complete history
         152             : !EOP
         153             : !------------------------------------------------------------------------------
         154             : !BOC
         155             : !
         156             : ! !DEFINED PARAMETERS:
         157             : !
         158             :     ! Universal gas constant for adjusting KOA for temp: 8.3144598 [J/mol/K]
         159             :     REAL(hp), PARAMETER :: R          = 8.3144598e+0_hp
         160             :
         161             :     ! Density of octanol, needed for partitioning into OC: 820 [kg/m^3]
         162             :     REAL(hp), PARAMETER :: DENS_OCT   = 82e+1_hp
         163             :
         164             :     ! Density of BC, needed for partitioning onto BC: 1 [kg/L] or 1000 [kg/m^3]
         165             :     ! From Lohmann and Lammel, Environ. Sci. Technol., 2004, 38:3793-3803.
         166             :     REAL(hp), PARAMETER :: DENS_BC    = 1e+3_hp
         167             : !
         168             : ! !LOCAL VARIABLES:
         169             : !
         170             :     INTEGER             :: I, J, L
         171             :     INTEGER             :: PBL_MAX
         172             :     INTEGER             :: MONTH,            YEAR
         173             :     REAL(hp)            :: F_OF_PBL,         TK
         174             :     REAL(hp)            :: T_POP
         175             :     REAL(hp)            :: C_OC1,            C_BC1
         176             :     REAL(hp)            :: C_OC2,            C_BC2
         177             :     REAL(hp)            :: F_POP_OC,         F_POP_BC
         178             :     REAL(hp)            :: F_POP_G,          AIR_VOL
         179             :     REAL(hp)            :: KOA_T,            KBC_T
         180             :     REAL(hp)            :: KOC_BC_T,         KBC_OC_T
         181             :     REAL(hp)            :: VR_OC_AIR,        VR_OC_BC
         182             :     REAL(hp)            :: VR_BC_AIR,        VR_BC_OC
         183             :     REAL(hp)            :: SUM_F
         184             :     REAL(hp)            :: OC_AIR_RATIO,     OC_BC_RATIO
         185             :     REAL(hp)            :: BC_AIR_RATIO,     BC_OC_RATIO
         186             :     REAL(hp)            :: FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND
         187             :     REAL(hp)            :: FRAC_LEAF, FRAC_LAKE, FRAC_SOIL
         188             : !    LOGICAL, SAVE       :: FIRST = .TRUE.
         189             :     LOGICAL             :: aIR
         190             :     LOGICAL             :: IS_SNOW_OR_ICE,   IS_LAND_OR_ICE
         191             :     CHARACTER(LEN=255)  :: MSG, LOC
         192             :
         193             :     ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
         194             :     ! from gas phase to OC. For now we use Delta H for phase transfer
         195             :     ! from the gas phase to the pure liquid state.
         196             :     ! For PHENANTHRENE:
         197             :     ! this is taken as the negative of the Delta H for phase transfer
         198             :     ! from the pure liquid state to the gas phase (Schwarzenbach,
         199             :     ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
         200             :     ! For PYRENE:
         201             :     ! this is taken as the negative of the Delta H for phase transfer
         202             :     ! from the pure liquid state to the gas phase (Schwarzenbach,
         203             :     ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
         204             :     ! For BENZO[a]PYRENE:
         205             :     ! this is also taken as the negative of the Delta H for phase transfer
         206             :     ! from the pure liquid state to the gas phase (Schwarzenbach,
         207             :     ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
         208             :     REAL(hp)            :: DEL_H
         209             :
         210             :     ! KOA_298 for partitioning of gas phase POP to atmospheric OC
         211             :     ! KOA_298 = Cpop in octanol/Cpop in atmosphere at 298 K
         212             :     ! For PHENANTHRENE:
         213             :     ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
         214             :     ! For PYRENE:
         215             :     ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
         216             :     ! For BENZO[a]PYRENE:
         217             :     ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
         218             :     ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
         219             :     REAL(hp)            :: KOA_298
         220             :
         221             :     ! KBC_298 for partitioning of gas phase POP to atmospheric BC
         222             :     ! KBC_298 = Cpop in black carbon/Cpop in atmosphere at 298 K
         223             :     ! For PHENANTHRENE:
         224             :     ! log KBC_298 = 10.0, or 1.0*10^10 [unitless]
         225             :     ! For PYRENE:
         226             :     ! log KBC_298 = 11.0, or 1.0*10^11 [unitless]
         227             :     ! For BENZO[a]PYRENE:
         228             :     ! log KBC_298 = 13.9, or 7.94*10^13 [unitless]
         229             :     ! (Lohmann and Lammel, EST, 2004, 38:3793-3802)
         230             :     REAL(hp)            :: KBC_298
         231             :
         232             :     ! Pointers
         233           0 :     REAL(sp),     POINTER :: Arr3D(:,:,:)
         234             :     TYPE(MyInst), POINTER :: Inst
         235             :
         236             :     !=======================================================================
         237             :     ! HCOX_GC_POPs_RUN begins here!
         238             :     !=======================================================================
         239           0 :     LOC = 'HCOX_GC_POPs_RUN (HCOX_GC_POPS_MOD.F90)'
         240             :
         241             :     ! Return if extension not turned on
         242           0 :     IF ( ExtState%GC_POPs <= 0 ) RETURN
         243             :
         244             :     ! Enter
         245           0 :     CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
         246           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         247           0 :         CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         248           0 :         RETURN
         249             :     ENDIF
         250             :
         251             :     ! Get instance
         252           0 :     Inst => NULL()
         253           0 :     CALL InstGet ( ExtState%GC_POPs, Inst, RC )
         254           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         255           0 :        WRITE(MSG,*) 'Cannot find GC_POPs instance Nr. ', ExtState%GC_POPs
         256           0 :        CALL HCO_ERROR(MSG,RC)
         257           0 :        RETURN
         258             :     ENDIF
         259             :
         260           0 :     DEL_H   = ExtState%POP_DEL_H
         261           0 :     KOA_298 = ExtState%POP_KOA
         262           0 :     KBC_298 = ExtState%POP_KBC
         263           0 :     Arr3D   => NULL()
         264             :
         265             :     !=======================================================================
         266             :     ! Get pointers to gridded data imported through config. file
         267             :     !=======================================================================
         268           0 :     IF ( HcoClock_First(HcoState%Clock,.TRUE.) ) THEN
         269             :
         270           0 :        CALL HCO_GetPtr( HcoState, 'TOT_POP',     Inst%POP_TOT_EM, RC )
         271           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         272           0 :            CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
         273           0 :            RETURN
         274             :        ENDIF
         275             :
         276           0 :        CALL HCO_GetPtr( HcoState, 'GLOBAL_OC',   Inst%C_OC,       RC )
         277           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         278           0 :            CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
         279           0 :            RETURN
         280             :        ENDIF
         281             :
         282           0 :        CALL HCO_GetPtr( HcoState, 'GLOBAL_BC',   Inst%C_BC,       RC )
         283           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         284           0 :            CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
         285           0 :            RETURN
         286             :        ENDIF
         287             :
         288           0 :        CALL HCO_GetPtr( HcoState, 'SURF_POP',    Inst%POP_SURF,   RC )
         289           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         290           0 :            CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
         291           0 :            RETURN
         292             :        ENDIF
         293             :
         294           0 :        CALL HCO_GetPtr( HcoState, 'SOIL_CARBON', Inst%F_OC_SOIL,  RC )
         295           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         296           0 :            CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
         297           0 :            RETURN
         298             :        ENDIF
         299             :
         300             :        ! Convert F_OC_SOIL from kg/m2 to fraction
         301           0 :        DO J=1, HcoState%NY
         302           0 :        DO I=1, HcoState%NX
         303             :
         304             :           ! Assume most of carbon mass extends to 5 cm and calculate
         305             :           ! concentration in kg/kg
         306             :           ! For now, assume a mean soil bulk density of 1300 kg/m3 similar to
         307             :           ! McLachlan 2002 to calculate a dry weight fraction
         308           0 :           Inst%F_OC_SOIL(I,J) = Inst%F_OC_SOIL(I,J) / 30e-2_hp / 13e+2_hp
         309             :
         310             :        ENDDO
         311             :        ENDDO
         312             :
         313             : !       FIRST = .FALSE.
         314             :
         315             :     ENDIF
         316             :
         317             :     ! Maximum extent of the PBL [model level]
         318           0 :     IF ( .NOT. ASSOCIATED(ExtState%PBL_MAX) ) THEN
         319           0 :        CALL HCO_ERROR ( 'PBL_MAX not defined in ExtState!', RC )
         320           0 :        RETURN
         321             :     ELSE
         322           0 :        PBL_MAX = DBLE( ExtState%PBL_MAX )
         323             :     ENDIF
         324             :
         325             :     !=================================================================
         326             :     ! Call emissions routines for revolatilization fluxes from surfaces
         327             :     ! Assume all re-emited POPs are in the gas phase until partitioning
         328             :     ! to ambient OC and BC in the boundary layer.
         329             :     ! Re-emission flux/mass depends on type of POP
         330             :     ! First draft, CLF, 28 Aug 2012
         331             :     !=================================================================
         332           0 :     CALL SOILEMISPOP( Inst%POP_SURF, Inst%F_OC_SOIL, HcoState, ExtState, Inst )
         333           0 :     CALL LAKEEMISPOP( Inst%POP_SURF, HcoState, ExtState, Inst )
         334           0 :     CALL VEGEMISPOP ( Inst%POP_SURF, HcoState, ExtState, Inst )
         335             :
         336           0 :     Inst%EmisPOPGfromSnow  = 0e+0_hp
         337           0 :     Inst%EmisPOPGfromOcean = 0e+0_hp
         338             :
         339             :     ! Loop over grid boxes
         340           0 :     DO J = 1, HcoState%Ny
         341           0 :     DO I = 1, HcoState%Nx
         342             :
         343           0 :        F_OF_PBL = 0e+0_hp
         344           0 :        T_POP    = 0e+0_hp
         345             :
         346             :        ! Here, save the total from the emissions array
         347             :        ! into the T_POP variable [kg/m2/s]
         348           0 :        T_POP = Inst%POP_TOT_EM(I,J)
         349             :
         350             :        ! Now add revolatilization (secondary) emissions to primary [kg/m2/s]
         351           0 :        T_POP = T_POP + Inst%EmisPOPGfromSoil(I,J) + &
         352           0 :                        Inst%EmisPOPGfromLake(I,J) + &
         353           0 :                        Inst%EmisPOPGfromLeaf(I,J) !+ &
         354             :                        !Inst%EmisPOPGfromSnow(I,J) + &
         355             :                        !Inst%EmisPOPGfromOcean(I,J)
         356             :
         357             :        !====================================================================
         358             :        ! Apportion total POPs emitted to gas phase, OC-bound, and BC-bound
         359             :        ! emissions (clf, 2/1/2011)
         360             :        ! Then partition POP throughout PBL; store into STT [kg]
         361             :        ! Now make sure STT does not underflow (cdh, bmy, 4/6/06; eck 9/20/10)
         362             :        !====================================================================
         363             :
         364             :        ! Loop up to max PBL level
         365           0 :        DO L = 1, PBL_MAX
         366             :
         367             :           ! Get temp [K]
         368           0 :           TK = ExtState%TK%Arr%Val(I,J,L)
         369             :
         370             :           ! Define temperature-dependent partition coefficients:
         371             :           ! KOA_T, the octanol-air partition coeff at temp T [unitless]
         372           0 :           KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp)))
         373             :
         374             :           ! Define KBC_T, the BC-air partition coeff at temp T [unitless]
         375             :           ! TURN OFF TEMPERATURE DEPENDENCY FOR SENSITIVITY ANALYSIS
         376           0 :           KBC_T = KBC_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp)))
         377             :
         378             :           ! Define KOC_BC_T, theoretical OC-BC part coeff at temp T [unitless]
         379           0 :           KOC_BC_T = KOA_T / KBC_T
         380             :
         381             :           ! Define KBC_OC_T, theoretical BC_OC part coeff at temp T [unitless]
         382           0 :           KBC_OC_T = 1d0 / KOC_BC_T
         383             :
         384             :           ! Get monthly mean OC and BC concentrations [kg/box]
         385           0 :           C_OC1    = Inst%C_OC(I,J,L)
         386           0 :           C_BC1    = Inst%C_BC(I,J,L)
         387             :
         388             :           ! Make sure OC is not negative
         389           0 :           C_OC1    = MAX( C_OC1, 0e+0_hp )
         390             :
         391             :           ! Convert C_OC and C_BC units to volume per box
         392             :           ! [m^3 OC or BC/box]
         393           0 :           C_OC2    = C_OC1 / DENS_OCT
         394           0 :           C_BC2    = C_BC1 / DENS_BC
         395             :
         396             :           ! Get air volume (m^3)
         397           0 :           AIR_VOL  = ExtState%AIRVOL%Arr%Val(I,J,L)
         398             :
         399             :           ! Define volume ratios:
         400             :           ! VR_OC_AIR = volume ratio of OC to air [unitless]
         401           0 :           VR_OC_AIR   = C_OC2 / AIR_VOL
         402             :
         403             :           ! VR_OC_BC  = volume ratio of OC to BC [unitless]
         404           0 :           VR_OC_BC    = C_OC2 / C_BC2
         405             :
         406             :           ! VR_BC_AIR = volume ratio of BC to air [unitless]
         407           0 :           VR_BC_AIR   = VR_OC_AIR / VR_OC_BC
         408             :
         409             :           ! VR_BC_OC  = volume ratio of BC to OC [unitless]
         410             :           !VR_BC_OC(I,J,L)    = 1d0 / VR_OC_BC(I,J,L)
         411           0 :           VR_BC_OC    = 1d0 / VR_OC_BC
         412             :
         413             :           ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase,
         414             :           ! BC-phase, and gas phase with new time step (should only change if
         415             :           ! temp changes or OC/BC concentrations change)
         416           0 :           OC_AIR_RATIO = 1e+0_hp / (KOA_T    * VR_OC_AIR)
         417           0 :           OC_BC_RATIO  = 1e+0_hp / (KOC_BC_T * VR_OC_BC)
         418             :
         419           0 :           BC_AIR_RATIO = 1e+0_hp / (KBC_T    * VR_BC_AIR)
         420           0 :           BC_OC_RATIO  = 1e+0_hp / (KBC_OC_T * VR_BC_OC)
         421             :
         422             :           ! If there are zeros in OC or BC concentrations, make sure they
         423             :           ! don't cause problems with phase fractions
         424           0 :           IF ( C_OC1 > SMALLNUM .and. C_BC1 > SMALLNUM ) THEN
         425           0 :              F_POP_OC  = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO + OC_BC_RATIO)
         426           0 :              F_POP_BC  = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO + BC_OC_RATIO)
         427             :
         428           0 :           ELSE IF ( C_OC1 > SMALLNUM .and. C_BC1 .le. SMALLNUM ) THEN
         429           0 :              F_POP_OC  = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO)
         430           0 :              F_POP_BC  = SMALLNUM
         431             :
         432           0 :           ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 > SMALLNUM ) THEN
         433           0 :              F_POP_OC  = SMALLNUM
         434           0 :              F_POP_BC  = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO)
         435             :
         436           0 :           ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 .le. SMALLNUM) THEN
         437           0 :              F_POP_OC = SMALLNUM
         438           0 :              F_POP_BC = SMALLNUM
         439             :           ENDIF
         440             :
         441             :           ! Gas-phase:
         442           0 :           F_POP_G   = 1e+0_hp - F_POP_OC - F_POP_BC
         443             :
         444             :           ! Check that sum of fractions equals 1
         445           0 :           SUM_F = F_POP_OC + F_POP_BC + F_POP_G
         446             :
         447             :           ! Fraction of PBL that box (I,J,L) makes up [unitless]
         448           0 :           F_OF_PBL = ExtState%FRAC_OF_PBL%Arr%Val(I,J,L)
         449             :
         450             :           ! Calculate rates of POP emissions in each phase [kg/m2/s]
         451             :           ! OC-phase:
         452           0 :           Inst%EmisPOPPOCPO(I,J,L) = F_POP_OC * F_OF_PBL * T_POP
         453             :
         454             :           ! BC-phase
         455           0 :           Inst%EmisPOPPBCPO(I,J,L) = F_POP_BC * F_OF_PBL * T_POP
         456             :
         457             :           ! Gas-phase
         458           0 :           Inst%EmisPOPG(I,J,L)  = F_POP_G  * F_OF_PBL * T_POP
         459             :
         460             :        ENDDO
         461             :
         462             :     ENDDO
         463             :     ENDDO
         464             :
         465             :     !=======================================================================
         466             :     ! Add POPs emissions to HEMCO data structure & diagnostics
         467             :     !=======================================================================
         468             :
         469             :     !----------------------
         470             :     ! OC-PHASE EMISSIONS
         471             :     !----------------------
         472           0 :     IF ( Inst%IDTPOPPOCPO > 0 ) THEN
         473             :
         474             :        ! Add flux to emissions array
         475           0 :        Arr3D => Inst%EmisPOPPOCPO(:,:,:)
         476             :        CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPOCPO, &
         477           0 :                          RC, ExtNr=Inst%ExtNr )
         478           0 :        Arr3D => NULL()
         479           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         480           0 :           CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPOCPO', RC )
         481           0 :           RETURN
         482             :        ENDIF
         483             :     ENDIF
         484             :
         485             :     !----------------------
         486             :     ! BC-PHASE EMISSIONS
         487             :     !----------------------
         488           0 :     IF ( Inst%IDTPOPPBCPO > 0 ) THEN
         489             :
         490             :        ! Add flux to emissions array
         491           0 :        Arr3D => Inst%EmisPOPPBCPO(:,:,:)
         492             :        CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPBCPO, &
         493           0 :                          RC, ExtNr=Inst%ExtNr )
         494           0 :        Arr3D => NULL()
         495           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         496           0 :           CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPBCPO', RC )
         497           0 :           RETURN
         498             :        ENDIF
         499             :     ENDIF
         500             :
         501             :     !----------------------
         502             :     ! GASEOUS EMISSIONS
         503             :     !----------------------
         504           0 :     IF ( Inst%IDTPOPG > 0 ) THEN
         505             :
         506             :        ! Add flux to emissions array
         507           0 :        Arr3D => Inst%EmisPOPG(:,:,:)
         508             :        CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPG, &
         509           0 :                          RC, ExtNr=Inst%ExtNr )
         510           0 :        Arr3D => NULL()
         511           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         512           0 :           CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPG', RC )
         513           0 :           RETURN
         514             :        ENDIF
         515             :
         516             :     ENDIF
         517             :
         518             :     !=======================================================================
         519             :     ! Cleanup & quit
         520             :     !=======================================================================
         521             :
         522             :     ! Nullify pointers
         523           0 :     Inst    => NULL()
         524             :
         525             :     ! Return w/ success
         526           0 :     CALL HCO_LEAVE( HcoState%Config%Err,RC )
         527             :
         528           0 :   END SUBROUTINE HCOX_GC_POPs_Run
         529             : !EOC
         530             : !------------------------------------------------------------------------------
         531             : !                  GEOS-Chem Global Chemical Transport Model                  !
         532             : !------------------------------------------------------------------------------
         533             : !BOP
         534             : !
         535             : ! !IROUTINE: soilemispop
         536             : !
         537             : ! !DESCRIPTION: Subroutine SOILEMISPOP is the subroutine for secondary
         538             : !  POP emissions from soils.
         539             : !\\
         540             : !\\
         541             : ! !INTERFACE:
         542             : !
         543           0 :       SUBROUTINE SOILEMISPOP( POP_SURF, F_OC_SOIL, HcoState, ExtState,  Inst )
         544             : !
         545             : ! !INPUT PARAMETERS:
         546             : !
         547             :       REAL(sp), DIMENSION(:,:), INTENT(IN)  :: POP_SURF   ! POP sfc conc [kg]
         548             :       REAL(sp), DIMENSION(:,:), INTENT(IN)  :: F_OC_SOIL  ! Frac C in soil [g/g]
         549             :       TYPE(HCO_STATE),          POINTER     :: HcoState   ! Hemco state
         550             :       TYPE(Ext_State),          POINTER     :: ExtState   ! Module options
         551             :       TYPE(MyInst),             POINTER     :: Inst       ! Instance
         552             : !
         553             : ! !REVISION HISTORY:
         554             : !  21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
         555             : !  See https://github.com/geoschem/hemco for complete history
         556             : !EOP
         557             : !------------------------------------------------------------------------------
         558             : !BOC
         559             : !
         560             : ! !LOCAL VARIABLES:
         561             : !
         562             :       LOGICAL  :: IS_SNOW_OR_ICE
         563             :       LOGICAL  :: IS_LAND_OR_ICE
         564             :       INTEGER  :: I, J, L
         565             :       REAL(hp) :: POPG
         566             :       REAL(hp) :: TK_SURF
         567             :       REAL(hp) :: KSA_T, FLUX, F_OC
         568             :       REAL(hp) :: SOIL_CONC, DTSRCE, KOA_T
         569             :       REAL(hp) :: DIFF, FSOIL, FAIR, DS
         570             :       REAL(hp) :: DSA, DAD, DWD, PL
         571             :       REAL(hp) :: ZSOIL, ZAIR, TK
         572             :       REAL(hp) :: DTCHEM, NEWSOIL, E_KDEG
         573             :       REAL(hp) :: FRAC_LAKE, FRAC_SOIL
         574             :       REAL(hp) :: FUG_R
         575             :       REAL(hp) :: AREA_M2
         576             : !
         577             : ! !DEFINED PARAMETERS:
         578             : !
         579             :       ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
         580             :       ! from gas phase to OC. For now we use Delta H for phase transfer
         581             :       ! from the gas phase to the pure liquid state.
         582             :       ! For PHENANTHRENE:
         583             :       ! this is taken as the negative of the Delta H for phase transfer
         584             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
         585             :       !  Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
         586             :       ! For PYRENE:
         587             :       ! this is taken as the negative of the Delta H for phase transfer
         588             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
         589             :       ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
         590             :       ! For BENZO[a]PYRENE:
         591             :       ! this is also taken as the negative of the Delta H for phase transfer
         592             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
         593             :       ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
         594             :       REAL(hp)            :: DEL_H
         595             :
         596             :       ! R = universal gas constant for adjusting KOA for temp:
         597             :       ! 8.3144598 [J/mol/K OR m3*Pa/K/mol]
         598             :       REAL(hp), PARAMETER :: R = 8.3144598d0
         599             :
         600             :       ! Molecular weight
         601             :       ! For phe, 0.17823 kg/mol
         602             :       ! For pyr, 0.20225 kg/mol
         603             :       ! For BaP, 0,25231 kg/mol
         604             :       REAL(hp)            :: MW
         605             :
         606             :       ! Molecular weight of air
         607             :       REAL(hp), PARAMETER :: MWAIR  = 28.97d0 ! g/mol
         608             :
         609             :       ! For PHENANTHRENE:
         610             :       ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
         611             :       ! For PYRENE:
         612             :       ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
         613             :       ! For BENZO[a]PYRENE:
         614             :       ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
         615             :       ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
         616             :       REAL(hp)            :: KOA_298
         617             :
         618             :       ! Set transfer velocity and diffusion coefficient values
         619             :       REAL(hp), PARAMETER :: KSA  = 1d0     ![m/h]
         620             :       REAL(hp), PARAMETER :: BA   = 0.04d0  ![m2/h]
         621             :       REAL(hp), PARAMETER :: BW   = 4d-6    ![m2/h]
         622             :
         623             :       ! Set soil degradation rate
         624             :       REAL(hp), PARAMETER :: DEGR = 3.5d-5  ![/h]
         625             :
         626             :       !=================================================================
         627             :       ! SOILEMISPOP begins here!
         628             :       !=================================================================
         629             :
         630             :       ! Copy values from ExtState
         631           0 :       DEL_H   = ExtState%POP_DEL_H
         632           0 :       KOA_298 = ExtState%POP_KOA
         633           0 :       MW      = ExtState%POP_XMW
         634             :
         635             :       ! Emission timestep [s]
         636           0 :       DTSRCE  = HcoState%TS_EMIS
         637             :
         638             :       ! Chemistry timestep [h]
         639           0 :       DTCHEM = HcoState%TS_CHEM / 60e+0_hp
         640             :
         641           0 :       DO J=1, HcoState%NY
         642           0 :       DO I=1, HcoState%NX
         643             :
         644             :          ! Set logicals
         645             :          ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
         646             :          ! IS_LAND will return non-ocean boxes but may still contain lakes
         647             :          ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
         648             :          IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR.  &
         649           0 :                             ( IS_ICE (I,J,ExtState) ) )
         650             :          IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR.  &
         651           0 :                             ( IS_LAND(I,J,ExtState)   .AND. &
         652           0 :                               ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
         653             :
         654             :          ! Do soils routine only if we are on land that is not covered with
         655             :          ! snow or ice
         656           0 :          IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
         657             :
         658             :             ! Get fraction of grid box covered by lake surface area
         659           0 :             FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J)
         660             :
         661             :             ! Get fraction of land remaining
         662             :             ! Assume the remaining land is soil and get OC content.
         663             :             ! If remaining land is not soil (e.g., desert), there
         664             :             ! should be a characteristically low OC content
         665             :             ! that will have little capacity to store POPs
         666             :             ! ONLY SUBTRACT FRAC LAKE NOW
         667           0 :             FRAC_SOIL = MAX(1e+0_hp - FRAC_LAKE, 0e+0_hp)
         668             :
         669             :             ! Get surface skin temp [K]
         670           0 :             TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
         671             :
         672             :             ! Get air temp [K]
         673           0 :             TK = ExtState%TK%Arr%Val(I,J,1)
         674             :
         675             :             ! Get gas phase air POP concentration at surface in mol/m3
         676             :             ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
         677           0 :             POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
         678             :
         679             :             ! (kg trc/kg dry air) / (0.178 kg trc/mol) * (kg dry air/m3)
         680           0 :             POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
         681             :
         682             :             ! Grid box surface area [m2]
         683           0 :             AREA_M2   = HcoState%Grid%AREA_M2%Val(I,J)
         684             :
         685             :             ! Get soil concentration in top 5 cm of soil
         686             :             ! (following Howsam et al 2000)
         687             :             ! From Howsam et al, soil burdens are equal to
         688             :             ! 2.6 years deposition for PHE,
         689             :             ! 10 years for PYR, and 9.4 years for BaP
         690             :             ! Convert to mol/m3
         691             :             ! 2.6 yrs * kg deposited to soil in 1 yr * / 0.178 kg/mol
         692             :             ! / area grid box (m2) / 0.05 m
         693           0 :             SOIL_CONC = 10e+0_hp * POP_SURF(I,J) / MW / &
         694           0 :                         AREA_M2 / 5e-2_hp ! mol/m3
         695             :
         696             :             ! Get rid of mass due to degradation
         697             :             ! Use rate constant for BaP from Mackay and Paterson 1991:
         698             :             ! 3.5*10^-5 /h
         699             :             ! Calculate exponential factor
         700           0 :             E_KDEG = EXP (-DEGR * DTCHEM)
         701             :
         702             :             ! Adjust conc
         703           0 :             NEWSOIL = SOIL_CONC * E_KDEG
         704             :
         705             :             ! Get foc from GTMM saved files
         706           0 :             F_OC = F_OC_SOIL(I,J)
         707             :
         708             :             ! Define temperature-dependent KOA:
         709             :             KOA_T = KOA_298 * EXP((-DEL_H/R) * ( ( 1e+0_hp / TK_SURF ) - &
         710           0 :                     ( 1e+0_hp / 298e+0_hp ) ))
         711             :
         712             :             ! Dimensionless coefficient (mol/m3 soil / mol/m3 air)
         713             :             ! KSA = 1.5 (fTOC)*Koa
         714           0 :             KSA_T = 1.5 * F_OC * KOA_T
         715           0 :             KSA_T = MAX( KSA_T, SMALLNUM )
         716             :
         717             :             ! Calculate fugacities from concentrations by dividing by "Z"
         718             :             ! values, or the fugacity capacity in mol/m3*Pa following
         719             :             ! Mackay and Paterson, 1991
         720             :
         721             :             ! fsoil = Csoil * R * T / KSA [Pa]
         722             :             ! where Csoil is in mol/m3, R is in Pa * m3 / mol K
         723             :             ! T is in K and KSA is dimensionless
         724           0 :             FSOIL = NEWSOIL * R * TK_SURF / KSA_T
         725             :
         726             :             ! fair = Cair * R * T [Pa]
         727             :             ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K
         728           0 :             FAIR = POPG * R * TK
         729             :
         730             :             ! Calculate the fugacity gradient [Pa]
         731             :             ! If the gradient is negative, fair is larger and the POP will
         732             :             ! diffuse from air to soil
         733             :             ! If the gradient is positive, fsoil is larger and POP will
         734             :             ! diffuse from soil to air
         735           0 :             DIFF  = FSOIL - FAIR
         736           0 :             FUG_R = FSOIL/FAIR
         737             :
         738             :             ! Calculate "Z" values from fugacities.
         739             :             ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f
         740           0 :             ZAIR  = POPG / FAIR ! (mol/m3) / (Pa)
         741           0 :             ZSOIL = NEWSOIL / FSOIL ! (mol/m3) / (Pa)
         742             :
         743             :             ! Calculate the "D" value, or the transfer coefficient that
         744             :             ! describes the movement of POP between phases (Mackay and
         745             :             ! Paterson, 1991). [mol/h*Pa]
         746             :             ! The D value for soil-air diffusion is given by
         747             :             ! Ds = 1 / (1/Dsa + 1/(Dad + Dwd))
         748             :             ! Dsa is the air-side boundary layer diffusion parameter [mol/h*Pa]
         749             :             ! Dad is the diffusion parameter between soil particles and
         750             :             !  "soil air" [mol/h*Pa]
         751             :             ! Dwd is the diffusion parameter between soil particles and
         752             :             !  porewater [mol/h*Pa]
         753             :             ! Dsa is in series with soil-air and soil-water diffusion,
         754             :             !  which are in parallel
         755             :
         756             :             ! Need to define each D value
         757             :             ! DSA = kSA * Zair
         758             :             ! where kSA is a mass transfer coefficient [m/h],
         759             :             ! Zair is the air fugacity capacity [mol/m3*Pa]
         760             :             ! ***********
         761             :             ! DAD = BA * Zair
         762             :             ! where BA is the molecular diffusivity in air [m2/h]
         763             :             ! ***********
         764             :             ! DWD = BW * Zsoil
         765             :             ! where BW is the molecular diffusivity in water [m2/h]
         766             :             ! **** PL = the soil diffusion pathlength, set to half the
         767             :             ! soil depth (0.025 m)
         768           0 :             DSA = KSA * ZAIR  ! (m/h)  * (mol/m3*Pa) = mol/m2*h*Pa
         769           0 :             DAD = BA* ZAIR    ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa
         770           0 :             DWD = BW * ZSOIL  ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa
         771           0 :             PL  = 0.025e+0_hp
         772           0 :             DS  = 1e+0_hp / ( 1e+0_hp/DSA + PL/(DAD+DWD) ) ! mol/(m2*h*Pa) [* m3*Pa/K/mol = m/h/K
         773             :
         774             :             ! Calculate Flux in mol/m2/h
         775           0 :             FLUX = DS * DIFF
         776             :
         777             :             ! Change to units of ng/m2/d for storage
         778           0 :             FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp
         779             :
         780             :             ! Kludge soil emissions from poles for now
         781             :             ! Bug somewhere that allows GCAP versions to think some high polar
         782             :             ! boxes during some months are land rather than ice - results in
         783             :             ! extremely high fluxes
         784           0 :             IF ( HcoState%Grid%YMID%Val(I,J) > 60    .OR. &
         785             :                  HcoState%Grid%YMID%Val(I,J) < -60 ) THEN
         786           0 :                FLUX = 0e+0_hp
         787           0 :                DIFF = 0e+0_hp
         788             :             ENDIF
         789             :
         790             :             ! Convert to an emission rate in kg/m2/s for returning to
         791             :             ! HcoX_GC_POPs_Run
         792           0 :             Inst%EmisPOPGfromSoil(I,J) = MAX( FLUX / 24e+0_hp / 3600e+0_hp / &
         793           0 :                                               1e+12_hp, 0e+0_hp )
         794             :
         795             :             ! Multiply the emissions by the fraction of land that is soil
         796           0 :             Inst%EmisPOPGfromSoil(I,J) = Inst%EmisPOPGfromSoil(I,J) * FRAC_SOIL
         797             :
         798             :             ! If the flux is positive, then the direction will be from the
         799             :             ! soil to the air.
         800             :             ! If the flux is zero or negative, store it in a separate array.
         801           0 :             IF ( FLUX > 0e+0_hp ) THEN
         802             :
         803             :                ! Store positive flux
         804           0 :                Inst%FluxPOPGfromSoilToAir(I,J) = FLUX
         805             :
         806             :                ! Make sure negative flux diagnostic has nothing added to it
         807           0 :                Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp
         808             :
         809             :                ! Store the soil/air fugacity ratio
         810           0 :                Inst%FugacitySoilToAir(I,J)   = FUG_R
         811             :
         812           0 :             ELSE IF ( FLUX <= 0e+0_hp ) THEN
         813             :
         814             :                ! Store the negative flux
         815           0 :                Inst%FluxPOPGfromAirToSoil(I,J) = FLUX
         816             :
         817             :                ! Add nothing to positive flux
         818           0 :                Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp
         819             :
         820             :                ! Continue to store the fugacity ratio
         821           0 :                Inst%FugacitySoilToAir(I,J)   = FUG_R
         822             :
         823             :             ENDIF
         824             :
         825             :          ELSE
         826             :
         827             :             ! We are not on land or the land is covered with ice or snow
         828           0 :             FLUX                            = 0e+0_hp
         829           0 :             Inst%EmisPOPGfromSoil(I,J)      = 0e+0_hp
         830           0 :             Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp
         831           0 :             Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp
         832           0 :             Inst%FugacitySoilToAir(I,J)      = 0e+0_hp
         833             :
         834             :          ENDIF
         835             :
         836             :       ENDDO
         837             :       ENDDO
         838             :
         839           0 :       END SUBROUTINE SOILEMISPOP
         840             : !EOC
         841             :
         842             : !------------------------------------------------------------------------------
         843             : !                  GEOS-Chem Global Chemical Transport Model                  !
         844             : !------------------------------------------------------------------------------
         845             : !BOP
         846             : !
         847             : ! !IROUTINE: lakeemispop
         848             : !
         849             : ! !DESCRIPTION: Subroutine LAKEEMISPOP is the subroutine for secondary
         850             : !  POP emissions from lakes.
         851             : !\\
         852             : !\\
         853             : ! !INTERFACE:
         854             : !
         855           0 :       SUBROUTINE LAKEEMISPOP( POP_SURF, HcoState, ExtState, Inst )
         856             : !
         857             : ! !INPUT PARAMETERS:
         858             : !
         859             :       REAL(sp), DIMENSION(:,:), INTENT(IN)  :: POP_SURF   ! POP sfc conc [kg]
         860             :       TYPE(HCO_STATE),          POINTER     :: HcoState   ! Hemco state
         861             :       TYPE(Ext_State),          POINTER     :: ExtState   ! Module options
         862             :       TYPE(MyInst),             POINTER     :: Inst       ! Instance
         863             : !
         864             : ! !REVISION HISTORY:
         865             : !  21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
         866             : !  See https://github.com/geoschem/hemco for complete history
         867             : !EOP
         868             : !------------------------------------------------------------------------------
         869             : !BOC
         870             : !
         871             : ! !LOCAL VARIABLES:
         872             : !
         873             :       INTEGER  :: I, J, L
         874             :       REAL(hp) :: TK_SURF
         875             :       REAL(hp) :: KAW_T, FLUX, KOL_T
         876             :       REAL(hp) :: DTSRCE
         877             :       REAL(hp) :: KA_H2O, KA_POP, KW_CO2, KW_POP
         878             :       REAL(hp) :: DA_H2O, DA_POP, DW_CO2, DW_POP
         879             :       REAL(hp) :: SCH_CO2, SCH_POP
         880             :       REAL(hp) :: TK, PRESS
         881             :       REAL(hp) :: C_DISS
         882             :       REAL(hp) :: POPG, U10M, ALPHA
         883             :       REAL(hp) :: FRAC_LAKE, VISC_H2O
         884             :       REAL(hp) :: SFCWINDSQR
         885             :       REAL(hp) :: AREA_M2
         886             :       LOGICAL  :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE
         887             : !
         888             : ! !DEFINED PARAMETERS:
         889             : !
         890             :       ! Delta H for POP:
         891             :       ! For PHENANTHRENE:
         892             :       ! this is the Delta H for phase transfer
         893             :       ! from air to water (Schwarzenbach,
         894             :       !  Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47 [kJ/mol].
         895             :       ! For PYRENE: 43000 [kJ/mol]
         896             :       REAL(hp)            :: DEL_HW
         897             :
         898             :       ! R = universal gas constant for adjusting KOA for temp:
         899             :       ! 8.3144598d-3 [kJ/mol/K]
         900             :       REAL(hp), PARAMETER :: R = 8.3144598d-3
         901             :
         902             :       ! Molecular weight
         903             :       ! For phe, 0.17823 kg/mol
         904             :       REAL(hp)            :: MWPOP
         905             :
         906             :       ! Molecular weight of air
         907             :       REAL(hp), PARAMETER :: MWAIR  = 28.97d0 ! g/mol
         908             :
         909             :       ! Molecular weight of water
         910             :       REAL(hp), PARAMETER :: MWH2O  = 18.1d0 ! g/mol
         911             :
         912             :       ! Molar volumes calculated following Abraham and McGowan 1987 as
         913             :       ! summarized by Schwarzenbach et al. 2003.
         914             :       ! Each element is assigned a characteristic atomic volume, and an atomic
         915             :       ! volume of 6.56 cm3/mol is subtracted for each bond, no matter whether
         916             :       ! single, double, or triple
         917             :       ! C = 16.35,  H = 8.71,  O = 12.43, N = 14.39, P = 24.87, F = 10.48
         918             :       ! Br = 26.21, I = 34.53, S = 22.91, Si = 26.83
         919             :
         920             :       ! Molar volume of water
         921             :       ! 2*(8.71) + 12.43 - 2*(6.56) = 16.73
         922             :       REAL(hp), PARAMETER :: V_H2O  = 16.73d0 ! cm3/mol
         923             :
         924             :       ! Molar volume of CO2
         925             :       ! 16.35 + 2*(12.43) - 2*(6.56) = 28.1
         926             :       REAL(hp), PARAMETER :: V_CO2   = 28.1d0  ! cm3/mol
         927             :
         928             :       ! Molar volume of POP
         929             :       ! For PHE (C16H10):
         930             :       ! 16*(16.35) + 10*(8.71) - 29*(6.56) =
         931             :       REAL(hp), PARAMETER :: V_POP  = 538.94d0 ! cm3/mol
         932             :
         933             :       ! Molar volume of air - average of gases in air
         934             :       REAL(hp), PARAMETER :: V_AIR  = 20.1d0     ! cm3/mol
         935             :
         936             :       ! For PHENANTHRENE:
         937             :       ! log KAW_298 = -2.76, or 1.74*10-3 [unitless]
         938             :       ! For PYRENE:
         939             :       ! log KAW_298 = -3.27, or 5.37*10-4 [unitless]
         940             :       REAL(hp)            :: KAW_298
         941             :
         942             :       ! Set the kinematic viscosity of freshwater at 20C
         943             :       !REAL(hp), PARAMETER :: VISC_H2O  ! = 1d0    ![cm2/s]
         944             :
         945             :       ! Set aqueous degradation rate
         946             :       !REAL(hp), PARAMETER :: DEGR      != 3.5d-5  ![/h]
         947             :
         948             :       !=================================================================
         949             :       ! LAKEEMISPOP begins here!
         950             :       !=================================================================
         951             :
         952             :       ! Copy values from ExtState
         953           0 :       DEL_HW  = ExtState%POP_DEL_HW
         954           0 :       MWPOP   = ExtState%POP_XMW
         955           0 :       KAW_298 = ExtState%POP_HSTAR
         956             :
         957             :       ! Emission timestep [s]
         958           0 :       DTSRCE  = HcoState%TS_EMIS
         959             :
         960           0 :       DO J=1, HcoState%NY
         961           0 :       DO I=1, HcoState%NX
         962             :
         963             :          ! Set logicals
         964             :          ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
         965             :          ! IS_LAND will return non-ocean boxes but may still contain lakes
         966             :          ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
         967             :          IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR.  &
         968           0 :                             ( IS_ICE (I,J,ExtState) ) )
         969             :          IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR.  &
         970           0 :                             ( IS_LAND(I,J,ExtState)   .AND. &
         971           0 :                               ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
         972             :
         973             :          ! Do soils routine only if we are on land that is not covered with
         974             :          ! snow or ice
         975           0 :          IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
         976             :
         977             :             ! Get fraction of grid box covered by lake surface area
         978           0 :             FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J)
         979             :
         980           0 :             IF ( FRAC_LAKE > 0e+0_hp ) THEN
         981             :
         982             :                ! Get surface skin temp [K]
         983           0 :                TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
         984             :
         985             :                ! Get air temp [K]
         986           0 :                TK = ExtState%TK%Arr%Val(I,J,1)
         987             :
         988             :                ! Get surface pressure at end of dynamic time step [hPa]
         989           0 :                PRESS = ExtState%PSC2_WET%Arr%Val(I,J)
         990             :
         991             :                ! Convert to units of atm
         992           0 :                PRESS = PRESS / 1013.25e+0_hp
         993             :
         994             :                ! Get gas phase air POP concentration at surface in mol/m3
         995             :                ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
         996           0 :                POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
         997             :
         998             :                ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3)
         999           0 :                POPG = POPG / MWPOP * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
        1000             :
        1001             :                ! Grid box surface area [m2]
        1002           0 :                AREA_M2   = HcoState%Grid%AREA_M2%Val(I,J)
        1003             :
        1004             :                ! Get the dissolved POP concentration at lake surface in mol/m3
        1005             :                ! Distribute the total deposited mass to a volume that best
        1006             :                ! matches observed dissolved concentrations
        1007             :                ! Future versions should consider aqueous particle concentrations
        1008             :                ! and sinking rates and photolytic/microbial degradation
        1009             :
        1010             :                ! Start with 1 m - scale by 100
        1011           0 :                C_DISS = POP_SURF(I,J) / MWPOP / AREA_M2 / 100e+0_hp
        1012             :
        1013             :                ! Wind speed at 10m altitude [m/s]
        1014           0 :                SFCWINDSQR = ExtState%U10M%Arr%Val(I,J)**2 &
        1015           0 :                           + ExtState%V10M%Arr%Val(I,J)**2
        1016           0 :                U10M       = SQRT( SFCWINDSQR )
        1017             :
        1018             :                ! Need to calculate water-side and air-side mass transfer
        1019             :                ! coefficients
        1020             :                ! Start with air-side
        1021             :                ! Relate air-side MTC of POP to that of H2O
        1022             :                ! First, calculate air-side MTC for water following
        1023             :                ! Schwarzenbach, Gschwend, Imboden 2003
        1024           0 :                KA_H2O = 0.2e+0_hp * U10M + 0.3e+0_hp   ! cm/s
        1025             :
        1026             :                ! Relate air-side MTC for water to that of POP via diffusivities
        1027             :                ! following Schwarzenbach et al 2003
        1028             :                ! Calculate temperature-dependent diffusivities in air first
        1029             :                ! folling Fuller et al. 1966 (summarized by Schwarzenbach
        1030             :                ! et al. 2003)
        1031             :                DA_POP = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + &
        1032             :                         (1e+0_hp / (MWPOP*1e+3_hp) ) )**0.5e+0_hp /    &
        1033             :                         ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) +         &
        1034           0 :                         V_POP**(1e+0_hp/3e+0_hp))**2e+0_hp )  ! cm2/s
        1035             :
        1036             :                DA_H2O = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + &
        1037             :                         (1e+0_hp / MWH2O ) )**0.5e+0_hp /              &
        1038             :                         ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) +         &
        1039           0 :                         V_H2O**(1e+0_hp/3e+0_hp))**2e+0_hp )  ! cm2/s
        1040             :
        1041             :                ! Relate POP and H2O air-side MTCs
        1042           0 :                KA_POP = KA_H2O * ( DA_POP / DA_H2O )**( 0.67e+0_hp ) ! cm/s
        1043             :
        1044             :                ! Now calculate water-side MTCs
        1045             :                ! Start with calculating the water side MTC of CO2
        1046             :                ! This depends on wind speed (Schwarzenbach et al 2003)
        1047             :                ! Three different scenarios for the water surface under
        1048             :                ! different wind speeds are considered:
        1049             :                ! Smooth Surface Regime (SSR): u10 <= 4.2 m/s
        1050             :                ! Rough Surface Regime (RSR):  4.2 m/s < u10 <= 13 m/s
        1051             :                ! Breaking Wave Regime (BWR):  u10 > 13 m/s
        1052             :                ! Alpha, the exponent in the relationship for the water side MTC,
        1053             :                ! also depends on the wind speed. Set this as well
        1054           0 :                IF ( U10M <= 4.2e+0_hp ) THEN
        1055             :                   KW_CO2 = 0.65e-3_hp       ! cm/s
        1056             :                   ALPHA  = 0.67e+0_hp
        1057           0 :                ELSE IF ( U10M > 4.2e+0_hp .AND. U10M <= 13e+0_hp ) THEN
        1058           0 :                   KW_CO2 = ( 0.79e+0_hp * U10M - 2.68e+0_hp ) * 1e-3_hp
        1059           0 :                   ALPHA  = 0.5e+0_hp
        1060           0 :                ELSE IF (U10M > 13e+0_hp) THEN
        1061           0 :                   KW_CO2 = ( 1.64e+0_hp * U10M - 13.69e+0_hp ) * 1e-3_hp
        1062           0 :                   ALPHA  = 0.50e+0_hp
        1063             :                ENDIF
        1064             :
        1065             :                ! Get the temperature-dependent kinematic viscosity of water
        1066           0 :                IF (TK_SURF <= 273.15 ) THEN
        1067             :                   VISC_H2O = 1.787e-2_hp      ! [cm2/s]
        1068           0 :                ELSE IF (TK_SURF > 273.15 .AND. TK_SURF <= 278.15 ) THEN
        1069             :                   VISC_H2O = 1.518e-2_hp
        1070           0 :                ELSE IF (TK_SURF > 258.15 .AND. TK_SURF <= 283.15 ) THEN
        1071             :                   VISC_H2O = 1.307e-2_hp
        1072           0 :                ELSE IF (TK_SURF > 283.15 .AND. TK_SURF <= 287.15 ) THEN
        1073             :                   VISC_H2O = 1.139e-2_hp
        1074           0 :                ELSE IF (TK_SURF > 287.15 .AND. TK_SURF <= 293.15 ) THEN
        1075             :                   VISC_H2O = 1.002e-2_hp
        1076           0 :                ELSE IF (TK_SURF > 293.15 .AND. TK_SURF <= 298.15 ) THEN
        1077             :                   VISC_H2O = 0.89e-2_hp
        1078           0 :                ELSE IF (TK_SURF > 298.15 ) THEN
        1079           0 :                   VISC_H2O = 0.797e-2_hp
        1080             :                ENDIF
        1081             :
        1082             :                ! Calculate the diffusivites of CO2 and POP in water
        1083             :                DW_CO2 = ( 13.26 * 1e-5_hp ) /               &
        1084             :                         (( VISC_H2O*1e+2_hp )**1.14e+0_hp * &
        1085           0 :                         (V_CO2)**0.589e+0_hp )   ! [cm2/s]
        1086             :
        1087             :                DW_POP = ( 13.26 * 1e-5_hp ) /               &
        1088             :                         (( VISC_H2O*1e+2_hp )**1.14e+0_hp * &
        1089           0 :                         (V_POP)**0.589e+0_hp )   ! [cm2/s]
        1090             :
        1091             :                ! Calculate the Schmidt numbers for CO2 and POP
        1092           0 :                SCH_CO2 = VISC_H2O / DW_CO2    ! [unitless]
        1093           0 :                SCH_POP = VISC_H2O / DW_POP    ! [unitless]
        1094             :
        1095             :                ! Calculate the water-side MTC for POP
        1096           0 :                KW_POP = KW_CO2 * ( SCH_POP / SCH_CO2 ) ** (-ALPHA) ! [cm/s]
        1097             :
        1098             :                ! Calculate the temperature-dependent dimensionless Henry's Law
        1099             :                ! constant
        1100             :                KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - &
        1101           0 :                         (1e+0_hp/298e+0_hp)))       ! [unitless]
        1102             :
        1103             :                ! Now calculate the overall air-water MTC
        1104             :                KOL_T = 1e+0_hp / ( 1e+0_hp/KW_POP + &
        1105           0 :                        1e+0_hp / (KA_POP*KAW_T) ) ! [cm/s]
        1106             :
        1107             :                ! Calculate Flux in ng/m2/day !
        1108             :                FLUX = KOL_T * 3600e+0_hp * 24e+0_hp * ( C_DISS - &
        1109           0 :                       POPG/KAW_T ) * MWPOP * 1e+12_hp / 100e+0_hp
        1110             :
        1111             :                ! Convert to an emission rate in kg/m2/s for returning to
        1112             :                ! HcoX_GC_POPs_Run
        1113             :                ! Only return it if it's positive
        1114           0 :                Inst%EmisPOPGfromLake(I,J) = MAX(FLUX / 24e+0_hp / 3600e+0_hp / &
        1115           0 :                                                 1e+12_hp, 0e+0_hp )
        1116             :
        1117             :                ! Multiply the emissions by the fraction of land that is water
        1118           0 :                Inst%EmisPOPGfromLake(I,J) = Inst%EmisPOPGfromLake(I,J) * &
        1119           0 :                                             FRAC_LAKE
        1120             :
        1121             :                ! If the flux is positive, then the direction will be from the
        1122             :                ! soil to the air.
        1123             :                ! If the flux is zero or negative, store it in a separate array.
        1124           0 :                IF ( FLUX > 0e+0_hp ) THEN
        1125             :
        1126             :                   ! Store positive flux
        1127           0 :                   Inst%FluxPOPGfromLakeToAir(I,J) = + FLUX
        1128             :
        1129             :                   ! Make sure negative flux diagnostic has nothing added to it
        1130           0 :                   Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp
        1131             :
        1132             :                      ! Store the soil/air fugacity ratio
        1133           0 :                   Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T)
        1134             :
        1135           0 :                ELSE IF ( FLUX <= 0e+0_hp ) THEN
        1136             :
        1137             :                   ! Store the negative flux
        1138           0 :                   Inst%FluxPOPGfromAirToLake(I,J) = FLUX
        1139             :
        1140             :                   ! Add nothing to positive flux
        1141           0 :                   Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp
        1142             :
        1143             :                   ! Continue to store the fugacity ratio
        1144           0 :                   Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T)
        1145             :
        1146             :                ENDIF
        1147             :
        1148             :             ENDIF
        1149             :
        1150             :          ELSE
        1151             :
        1152             :             ! We are not on land or the land is covered with ice or snow
        1153             :             ! or we are land but there is no water
        1154           0 :             FLUX                            = 0e+0_hp
        1155           0 :             Inst%EmisPOPGfromLake           = 0e+0_hp
        1156           0 :             Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp
        1157           0 :             Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp
        1158           0 :             Inst%FugacityLakeToAir(I,J)     = 0e+0_hp
        1159             :
        1160             :          ENDIF
        1161             :
        1162             :       ENDDO
        1163             :       ENDDO
        1164             :
        1165           0 :       END SUBROUTINE LAKEEMISPOP
        1166             : !EOC
        1167             : !------------------------------------------------------------------------------
        1168             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1169             : !------------------------------------------------------------------------------
        1170             : !BOP
        1171             : !
        1172             : ! !IROUTINE: vegemispop
        1173             : !
        1174             : ! !DESCRIPTION: Subroutine VEGEMISPOP is the subroutine for secondary
        1175             : !  POP emissions from leaves.
        1176             : !\\
        1177             : !\\
        1178             : ! !INTERFACE:
        1179             : !
        1180           0 :       SUBROUTINE VEGEMISPOP( POP_SURF, HcoState, ExtState, Inst )
        1181             : !
        1182             : ! !INPUT PARAMETERS:
        1183             : !
        1184             :       REAL(sp), DIMENSION(:,:), INTENT(IN)  :: POP_SURF   ! POP sfc conc [kg]
        1185             :       TYPE(HCO_State),          POINTER     :: HcoState   ! Hemco state
        1186             :       TYPE(Ext_State),          POINTER     :: ExtState   ! Module options
        1187             :       TYPE(MyInst),             POINTER     :: Inst       ! Instance
        1188             : !
        1189             : ! !REVISION HISTORY:
        1190             : !  21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
        1191             : !  See https://github.com/geoschem/hemco for complete history
        1192             : !EOP
        1193             : !------------------------------------------------------------------------------
        1194             : !BOC
        1195             : !
        1196             : ! !LOCAL VARIABLES:
        1197             : !
        1198             :       INTEGER  :: I, J, L
        1199             :       REAL(hp) :: POPG, POPG_GL
        1200             :       REAL(hp) :: FRAC_SNOWFREE_LAND, TK_SURF
        1201             :       REAL(hp) :: KLA_T, FLUX, K_MT
        1202             :       REAL(hp) :: LEAF_CONC, DTSRCE, KOA_T
        1203             :       REAL(hp) :: DIFF, FLEAF, FAIR, DS
        1204             :       REAL(hp) :: DAB_F, DC, DAL, PC, UC
        1205             :       REAL(hp) :: ZLEAF, ZAIR, TK
        1206             :       REAL(hp) :: NEWLEAF
        1207             :       REAL(hp) :: LAI, KAW_T, KOW_T
        1208             :       REAL(hp) :: FUG_R, NEW_LEAF, DLA
        1209             :       REAL(hp) :: AREA_M2
        1210             :       LOGICAL  :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE
        1211             :       LOGICAL  :: IS_SNOWFREE_LAND
        1212             :
        1213             : !
        1214             : ! !DEFINED PARAMETERS:
        1215             : !
        1216             :       ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
        1217             :       ! from gas phase to OC. For now we use Delta H for phase transfer
        1218             :       ! from the gas phase to the pure liquid state.
        1219             :       ! For PHENANTHRENE:
        1220             :       ! this is taken as the negative of the Delta H for phase transfer
        1221             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
        1222             :       !  Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
        1223             :       ! For PYRENE:
        1224             :       ! this is taken as the negative of the Delta H for phase transfer
        1225             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
        1226             :       ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
        1227             :       ! For BENZO[a]PYRENE:
        1228             :       ! this is also taken as the negative of the Delta H for phase transfer
        1229             :       ! from the pure liquid state to the gas phase (Schwarzenbach,
        1230             :       ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
        1231             :       REAL(hp)             :: DEL_H
        1232             :
        1233             :       ! Delta H for POP [kJ/mol].
        1234             :       ! For PHENANTHRENE:
        1235             :       ! this is the Delta H for phase transfer
        1236             :       ! from air to water (Schwarzenbach,
        1237             :       !  Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47000 [J/mol].
        1238             :       ! For PYRENE: 43000 [J/mol]
        1239             :       REAL(hp)            :: DEL_HW
        1240             :
        1241             :       ! R = universal gas constant for adjusting KOA for temp:
        1242             :       ! 8.314 [J/mol/K OR m3*Pa/K/mol]
        1243             :       REAL(hp), PARAMETER :: R = 8.3144598e+0_hp
        1244             :
        1245             :       ! Molecular weight
        1246             :       ! For phe, 0.17823 kg/mol
        1247             :       REAL(hp)            :: MW
        1248             :
        1249             :       ! Molecular weight of air
        1250             :       REAL(hp), PARAMETER :: MWAIR  = 28.97d0 ! g/mol
        1251             :
        1252             :       ! For PHENANTHRENE:
        1253             :       ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
        1254             :       ! For PYRENE:
        1255             :       ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
        1256             :       ! For BENZO[a]PYRENE:
        1257             :       ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
        1258             :       ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
        1259             :       REAL(hp)            :: KOA_298
        1260             :
        1261             :       ! For PHENANTHRENE:
        1262             :       ! log KAW_298 = -2.76, or 1.74*10-3 [unitless]
        1263             :       ! For PYRENE:
        1264             :       ! log KAW_298 = -3.27, or 5.37*10-4 [unitless]
        1265             :       REAL(hp)            :: KAW_298
        1266             :
        1267             :       ! Set volume fractions of octanol and water in surface and reservoir
        1268             :       ! leaf compartments [unitless]
        1269             :       REAL(hp), PARAMETER :: OCT_SURF = 0.8e+0_hp
        1270             :       REAL(hp), PARAMETER :: OCT_RES  = 0.02e+0_hp
        1271             :       REAL(hp), PARAMETER :: H2O_RES  = 0.7e+0_hp
        1272             :
        1273             :       ! Set thickness of different leaf compartments. Volumes calculated by
        1274             :       ! multiplying thicknesses by leaf area index
        1275             :       REAL(hp), PARAMETER :: SURF_THICK = 2e-6_hp   ! m
        1276             :       REAL(hp), PARAMETER :: RES_THICK  = 250e-6_hp ! m
        1277             :
        1278             :       ! Set transfer velocity and diffusion coefficient values
        1279             :       REAL(hp), PARAMETER :: UAB_F = 9e+0_hp   ![m/h]
        1280             :
        1281             :       ! Set soil degradation rate
        1282             :       REAL(hp), PARAMETER :: DEGR  = 3.5e-5_hp ![/h]
        1283             :
        1284             :       !=================================================================
        1285             :       ! VEGEMISPOP begins here!
        1286             :       !=================================================================
        1287             :
        1288             :       ! Copy values from ExtState
        1289           0 :       DEL_H   = ExtState%POP_DEL_H
        1290           0 :       KOA_298 = ExtState%POP_KOA
        1291           0 :       DEL_HW  = ExtState%POP_DEL_Hw
        1292           0 :       KAW_298 = ExtState%POP_HSTAR
        1293           0 :       MW      = ExtState%POP_XMW
        1294             :
        1295             :       ! Emission timestep [s]
        1296           0 :       DTSRCE  = HcoState%TS_EMIS
        1297             :
        1298           0 :       DO J=1, HcoState%NY
        1299           0 :       DO I=1, HcoState%NX
        1300             :
        1301             :          ! Set logicals
        1302             :          ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
        1303             :          ! IS_LAND will return non-ocean boxes but may still contain lakes
        1304             :          ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
        1305             :          IS_LAND_OR_ICE = ( (IS_LAND(I,J,ExtState)) .OR.  &
        1306           0 :                             (IS_ICE (I,J,ExtState)) )
        1307             :          IS_SNOW_OR_ICE = ( (IS_ICE (I,J,ExtState)) .OR.  &
        1308           0 :                             (IS_LAND(I,J,ExtState)  .AND. &
        1309           0 :                              ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
        1310             :
        1311             :          ! Do soils routine only if we are on land that is not covered with
        1312             :          ! snow or ice
        1313           0 :          IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
        1314             :
        1315             :             ! Get fraction of grid box covered by leaf surface area
        1316             :             ! Do not consider different vegetation types for now
        1317           0 :             LAI = ExtState%LAI%Arr%Val(I,J)
        1318             :
        1319           0 :             IF ( LAI > 0e+0_hp ) THEN
        1320             :
        1321             :                ! Get surface skin temp [K]
        1322           0 :                TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
        1323             :
        1324             :                ! Get air temp [K]
        1325           0 :                TK = ExtState%TK%Arr%Val(I,J,1)
        1326             :
        1327             :                ! Get gas phase air POP concentration at surface in mol/m3
        1328             :                ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
        1329           0 :                POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
        1330             :
        1331             :                ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3)
        1332           0 :                POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
        1333             :
        1334             :                ! Grid box surface area [m2]
        1335           0 :                AREA_M2   = HcoState%Grid%AREA_M2%Val(I,J)
        1336             :
        1337             :                ! Only consider partitioning into leaf surface (not reservoir)
        1338             :                ! for now following Mackay et al 2006 Environ Sci & Pollut Res
        1339             :                ! Include reservoir when land-atm models become dynamic
        1340             :
        1341             :                ! Assume that all leaf surfaces contain an average lipid content
        1342             :                ! of 80% (Mackay et al 2006)
        1343             :
        1344             :                ! Get leaf concentration
        1345             :                ! Convert to mol/m3
        1346             :                ! kg deposited to leaf in 1 yr * / 0.178 kg/mol
        1347             :                ! / area grid box (m2) / surface thickness m
        1348           0 :                LEAF_CONC = POP_SURF(I,J) / MW / AREA_M2 / SURF_THICK ! mol/m3
        1349             :
        1350             :                ! Check concentration in leaves by assuming a density similar
        1351             :                ! to water (1 g/cm3)
        1352             :
        1353             :                ! No degradation/metabolism for now. Just scale leaf
        1354             :                ! concentrations to match flux observations
        1355           0 :                NEWLEAF = LEAF_CONC/1e+4_hp  !SCALING FACTOR
        1356             :
        1357             :                ! Define temperature-dependent KOA:
        1358             :                KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK_SURF) - &
        1359           0 :                        (1e+0_hp/298e+0_hp)))
        1360             :
        1361             :                ! Calculate the temperature-dependent dimensionless Henry's Law
        1362             :                ! constant
        1363             :                KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - &
        1364           0 :                        (1e+0_hp/298e+0_hp)))       ! [unitless]
        1365             :
        1366             :                ! Estimate the temperature-dependent dimensionless octanol-water
        1367             :                ! constant
        1368           0 :                KOW_T = KOA_T * KAW_T
        1369             :
        1370             :                ! Define dimensionless leaf surface-air partition coefficient
        1371             :                ! (mol/m3 leaf / mol/m3 air)
        1372             :                ! KLA = foct_surf * Koa
        1373           0 :                KLA_T = OCT_SURF * KOA_T
        1374             : !               KLA_T = MAX( KLA_T, SMALLNUM )
        1375             :
        1376             :                ! Calculate fugacities from concentrations by dividing by "Z"
        1377             :                ! values, or the fugacity capacity in mol/m3*Pa following Mackay
        1378             :                ! and Paterson, 1991
        1379             :
        1380             :                ! fleaf = Cleaf * R * T / KSA [Pa]
        1381             :                ! where Cleaf is in mol/m3, R is in Pa * m3 / mol K
        1382             :                ! T is in K and KLA is dimensionless
        1383           0 :                FLEAF = NEWLEAF * R * TK_SURF / KLA_T
        1384             :
        1385             :                ! fair = Cair * R * T [Pa]
        1386             :                ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K
        1387           0 :                FAIR = POPG * R * TK
        1388             :
        1389             :                ! Calculate the fugacity gradient [Pa]
        1390             :                ! If the gradient is negative, fair is larger and the POP will
        1391             :                ! diffuse from air to soil
        1392             :                ! If the gradient is positive, fsoil is larger and POP will
        1393             :                ! diffuse from soil to air
        1394           0 :                DIFF = FLEAF - FAIR
        1395           0 :                FUG_R = FLEAF/FAIR
        1396             :
        1397             :                ! Calculate "Z" values from fugacities.
        1398             :                ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f
        1399           0 :                ZAIR  = POPG / FAIR ! (mol/m3) / (Pa)
        1400           0 :                ZLEAF = NEWLEAF / FLEAF ! (mol/m3) / (Pa)
        1401             :
        1402             :                ! Calculate the "D" value, or the transfer coefficient that
        1403             :                ! describes the movement of POP between phases (Mackay and
        1404             :                ! Paterson, 1991, Cousins and Mackay 2000, internal report).
        1405             :                ! [mol/h*Pa]
        1406             :                ! The D value for leaf surface-air gas diffusion is given by
        1407             :                ! Dla = 1 / (1/Dc + 1/(Dab-f)) [mol/(Pa*h)]
        1408             :                ! where Dab-f is the boundary layer diffusion  [mol/h*Pa]
        1409             :                ! given by Dab-f = As * L * Uab-f * Za
        1410             :                ! where As is the area of the land surface [m2], L is the leaf
        1411             :                ! area index [m2/m2],
        1412             :                ! Uab-f is a mass transfer coefficient for surface-air boundary
        1413             :                ! layer diffusion [m/h],
        1414             :                ! and Za is the fugacity capacity of the air [mol/(m3*Pa)]
        1415             :                ! Dc is the cuticle diffusion, given by
        1416             :                ! Dc = As * L * Uc * Zf
        1417             :                ! where As and L are as above, Uc is the cuticle mass transfer
        1418             :                ! coefficient [m/h],
        1419             :                ! and Zf is the fugacity capacity of the leaf surface
        1420             :                ! (mol/(m3*Pa))
        1421             :
        1422             :                ! Uc is given by
        1423             :                ! Uc = 3600 * Pc * 1/Kaw
        1424             :                ! where Pc is the cuticle permeance (m/s) and Kaw is the
        1425             :                ! dimensionless air-water partition coefficient.
        1426             :                ! Pc is given by
        1427             :                ! Log Pc = ((0.704 * log Kow - 11.2) +
        1428             :                !          (-3.47 - 2.79 * logMW + 0.970 log Kow)) / 2
        1429             :                ! (an average of two equations)
        1430             :
        1431             :                ! Need to define each D value
        1432             :                ! DAB_F:
        1433             :                !  m/h * mol/(m3*Pa)  =  (mol/h*Pa*m2)
        1434           0 :                DAB_F = UAB_F * ZAIR  !  mol/(h*Pa*m2)
        1435             :
        1436             :                ! Calculate PC and then Uc in order to calculate Dc
        1437             :                ! PC, UC are calculated according to Cousins and Mackay,
        1438             :                ! Chemosphere, 2001, Table 2
        1439             :                PC = 10** (( 0.704e+0_hp * LOG (KOW_T) - 11.2e+0_hp ) +    &
        1440             :                     ( -3.47e+0_hp -2.79e+0_hp* LOG(MW*1000d0) + 0.97e+0_hp &
        1441           0 :                     * LOG(KOW_T))/2e+0_hp) ![m/s]
        1442             :
        1443           0 :                UC = 3600e+0_hp * PC * 1e+0_hp/KAW_T ! [m/h]
        1444             :
        1445           0 :                DC = UC * ZLEAF    ! mol/(h*Pa*m2)
        1446             :
        1447             :                ! Now calculate overall transfer  ! mol/(h*Pa*m2)
        1448           0 :                DLA = 1e+0_hp / (1e+0_hp/DC + 1e+0_hp/DAB_F)
        1449             :
        1450             :                ! Calculate Flux in mol/h/m2
        1451           0 :                FLUX = DLA * DIFF
        1452             :
        1453             :                ! Change to units of ng/m2/d for storage
        1454           0 :                FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp
        1455             :
        1456             :                ! Convert to an emission rate in kg/m2/s for returning to
        1457             :                ! HcoX_GC_POPs_Run
        1458             :                ! Only want to add rates that are positive
        1459           0 :                Inst%EmisPOPGfromLeaf(I,J) = MAX(FLUX * LAI / 24e+0_hp / &
        1460           0 :                                                 3600e+0_hp / 1e+12_hp, 0e+0_hp)
        1461             :
        1462             :                ! If the flux is positive, then the direction will be from the
        1463             :                ! soil to the air.
        1464             :                ! If the flux is zero or negative, store it in a separate array.
        1465           0 :                IF ( FLUX > 0e+0_hp ) THEN
        1466             :
        1467             :                   ! Store positive flux
        1468           0 :                   Inst%FluxPOPGfromLeafToAir(I,J) = FLUX
        1469             :
        1470             :                   ! Make sure negative flux diagnostic has nothing added to it
        1471           0 :                   Inst%FluxPOPGfromAirtoLeaf = 0e+0_hp
        1472             :
        1473             :                   ! Store the soil/air fugacity ratio
        1474           0 :                   Inst%FugacityLeafToAir(I,J) = FUG_R
        1475             :
        1476           0 :                ELSE IF ( FLUX <= 0e+0_hp ) THEN
        1477             :
        1478             :                   ! Store the negative flux
        1479           0 :                   Inst%FluxPOPGfromAirtoLeaf(I,J) = FLUX
        1480             :
        1481             :                   ! Add nothing to positive flux
        1482           0 :                   Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp
        1483             :
        1484             :                   ! Continue to store the fugacity ratio
        1485           0 :                   Inst%FugacityLeafToAir(I,J)   = FUG_R
        1486             :
        1487             :                ENDIF
        1488             :
        1489             :             ELSE
        1490             :
        1491             :                ! We are not on land or the land is covered with ice or snow
        1492           0 :                FLUX                            = 0e+0_hp
        1493           0 :                Inst%EmisPOPGfromLeaf(I,J)      = 0e+0_hp
        1494           0 :                Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp
        1495           0 :                Inst%FluxPOPGfromAirtoLeaf(I,J) = 0e+0_hp
        1496           0 :                Inst%FugacityLeafToAir(I,J)     = 0e+0_hp
        1497             :
        1498             :             ENDIF
        1499             :
        1500             :          ENDIF
        1501             :
        1502             :       ENDDO
        1503             :       ENDDO
        1504             :
        1505           0 :       END SUBROUTINE VEGEMISPOP
        1506             : !EOC
        1507             : !------------------------------------------------------------------------------
        1508             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1509             : !------------------------------------------------------------------------------
        1510             : !BOP
        1511             : !
        1512             : ! !IROUTINE: is_land
        1513             : !
        1514             : ! !DESCRIPTION: Function IS\_LAND returns TRUE if surface grid box (I,J) is
        1515             : !  a land box.
        1516             : !\\
        1517             : !\\
        1518             : ! !INTERFACE:
        1519             : !
        1520           0 :       FUNCTION IS_LAND( I, J, ExtState ) RESULT ( LAND )
        1521             : !
        1522             : ! !USES:
        1523             : !
        1524             :     USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE
        1525             : !
        1526             : ! !INPUT PARAMETERS:
        1527             : !
        1528             :       INTEGER,         INTENT(IN) :: I           ! Longitude index of grid box
        1529             :       INTEGER,         INTENT(IN) :: J           ! Latitude  index of grid box
        1530             :       TYPE(Ext_State), POINTER    :: ExtState    ! Module options
        1531             : !
        1532             : ! !RETURN VALUE:
        1533             : !
        1534             :       LOGICAL                     :: LAND        ! =T if it is a land box
        1535             : !
        1536             : ! !REVISION HISTORY:
        1537             : !  26 Jun 2000 - R. Yantosca - Initial version
        1538             : !  See https://github.com/geoschem/hemco for complete history
        1539             : !EOP
        1540             : !------------------------------------------------------------------------------
        1541             : !BOC
        1542           0 :       LAND = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J),   &
        1543           0 :                            ExtState%FRLANDIC%Arr%Val(I,J), &
        1544           0 :                            ExtState%FROCEAN%Arr%Val(I,J),  &
        1545           0 :                            ExtState%FRSEAICE%Arr%Val(I,J), &
        1546           0 :                            ExtState%FRLAKE%Arr%Val(I,J) ) == 1
        1547             :
        1548           0 :       END FUNCTION IS_LAND
        1549             : !EOC
        1550             : !------------------------------------------------------------------------------
        1551             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1552             : !------------------------------------------------------------------------------
        1553             : !BOP
        1554             : !
        1555             : ! !IROUTINE: is_ice
        1556             : !
        1557             : ! !DESCRIPTION: Function IS\_ICE returns TRUE if surface grid box (I,J)
        1558             : !  contains either land-ice or sea-ice.
        1559             : !\\
        1560             : !\\
        1561             : ! !INTERFACE:
        1562             : !
        1563           0 :       FUNCTION IS_ICE( I, J, ExtState ) RESULT ( ICE )
        1564             : !
        1565             : ! !USES:
        1566             : !
        1567             :     USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE
        1568             : !
        1569             : ! !INPUT PARAMETERS:
        1570             : !
        1571             :       INTEGER,         INTENT(IN) :: I           ! Longitude index of grid box
        1572             :       INTEGER,         INTENT(IN) :: J           ! Latitude  index of grid box
        1573             :       TYPE(Ext_State), POINTER    :: ExtState    ! Module options
        1574             : !
        1575             : ! !RETURN VALUE:
        1576             : !
        1577             :       LOGICAL                     :: ICE         ! =T if this is an ice box
        1578             : !
        1579             : ! !REVISION HISTORY:
        1580             : !  09 Aug 2005 - R. Yantosca - Initial version
        1581             : !  See https://github.com/geoschem/hemco for complete history
        1582             : !EOP
        1583             : !------------------------------------------------------------------------------
        1584             : !BOC
        1585             :
        1586           0 :       ICE = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J),   &
        1587           0 :                           ExtState%FRLANDIC%Arr%Val(I,J), &
        1588           0 :                           ExtState%FROCEAN%Arr%Val(I,J),  &
        1589           0 :                           ExtState%FRSEAICE%Arr%Val(I,J), &
        1590           0 :                           ExtState%FRLAKE%Arr%Val(I,J) ) == 2
        1591             :
        1592           0 :       END FUNCTION IS_ICE
        1593             : !EOC
        1594             : !------------------------------------------------------------------------------
        1595             : !                   Harmonized Emissions Component (HEMCO)                    !
        1596             : !------------------------------------------------------------------------------
        1597             : !BOP
        1598             : !
        1599             : ! !IROUTINE: HCOX_GC_POPs_Init
        1600             : !
        1601             : ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Init initializes the HEMCO
        1602             : ! GC\_POPs extension.
        1603             : !\\
        1604             : !\\
        1605             : ! !INTERFACE:
        1606             : !
        1607           0 :   SUBROUTINE HCOX_GC_POPs_Init( HcoState, ExtName, ExtState, RC )
        1608             : !
        1609             : ! !USES:
        1610             : !
        1611             :     USE HCO_ExtList_Mod,   ONLY : GetExtNr
        1612             :     USE HCO_STATE_MOD,     ONLY : HCO_GetExtHcoID
        1613             : !
        1614             : ! !INPUT PARAMETERS:
        1615             : !
        1616             :     CHARACTER(LEN=*), INTENT(IN   )  :: ExtName     ! Extension name
        1617             :     TYPE(Ext_State),  POINTER        :: ExtState    ! Module options
        1618             :     TYPE(HCO_State),  POINTER        :: HcoState    ! Hemco state
        1619             : !
        1620             : ! !INPUT/OUTPUT PARAMETERS:
        1621             : !
        1622             :     INTEGER,          INTENT(INOUT)  :: RC
        1623             :
        1624             : ! !REVISION HISTORY:
        1625             : !  19 Aug 2014 - M. Sulprizio- Initial version
        1626             : !  See https://github.com/geoschem/hemco for complete history
        1627             : !EOP
        1628             : !------------------------------------------------------------------------------
        1629             : !BOC
        1630             : !
        1631             : ! !LOCAL VARIABLES:
        1632             : !
        1633             :     ! Scalars
        1634             :     INTEGER                        :: I, N, nSpc, ExtNr
        1635             :     CHARACTER(LEN=255)             :: MSG, LOC
        1636             :     CHARACTER(LEN=255)             :: DiagnName
        1637             :     CHARACTER(LEN=255)             :: OutUnit
        1638             :
        1639             :     ! Arrays
        1640           0 :     INTEGER,           ALLOCATABLE :: HcoIDs(:)
        1641           0 :     CHARACTER(LEN=31), ALLOCATABLE :: SpcNames(:)
        1642             :
        1643             :     ! Pointers
        1644             :     TYPE(MyInst), POINTER          :: Inst
        1645           0 :     REAL(sp),     POINTER          :: Target2D(:,:)
        1646             :
        1647             :     !=======================================================================
        1648             :     ! HCOX_GC_POPs_INIT begins here!
        1649             :     !=======================================================================
        1650           0 :     LOC = 'HCOX_GC_POPs_INIT (HCOX_GC_POPS_MOD.F90)'
        1651             :
        1652             :     ! Get the extension number
        1653           0 :     ExtNr = GetExtNr( HcoState%Config%ExtList, TRIM(ExtName) )
        1654           0 :     IF ( ExtNr <= 0 ) RETURN
        1655             :
        1656             :     ! Enter HEMCO
        1657           0 :     CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
        1658           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1659           0 :         CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
        1660           0 :         RETURN
        1661             :     ENDIF
        1662             :
        1663             :     ! Create Instance
        1664           0 :     Inst => NULL()
        1665           0 :     CALL InstCreate ( ExtNr, ExtState%GC_POPs, Inst, RC )
        1666           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1667           0 :        CALL HCO_ERROR ( 'Cannot create GC_POPs instance', RC )
        1668           0 :        RETURN
        1669             :     ENDIF
        1670             :     ! Also fill ExtNrSS - this is the same as the parent ExtNr
        1671           0 :     Inst%ExtNr = ExtNr
        1672             :
        1673             :     ! Set species IDs
        1674           0 :     CALL HCO_GetExtHcoID( HcoState, Inst%ExtNr, HcoIDs, SpcNames, nSpc, RC )
        1675           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1676           0 :         CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
        1677           0 :         RETURN
        1678             :     ENDIF
        1679             :
        1680             :     ! Verbose mode
        1681           0 :     IF ( HcoState%amIRoot ) THEN
        1682             :
        1683             :        ! Write the name of the extension regardless of the verbose setting
        1684           0 :        msg = 'Using HEMCO extension: GC_POPs (POPs emissions)'
        1685           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        1686           0 :           CALL HCO_Msg( HcoState%Config%Err, msg, sep1='-' ) ! with separator
        1687             :        ELSE
        1688           0 :           CALL HCO_Msg( msg, verb=.TRUE.                   ) ! w/o separator
        1689             :        ENDIF
        1690             :
        1691             :        ! Write all other messages as debug printout only
        1692           0 :        MSG = 'Use the following species (Name: HcoID):'
        1693           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        1694           0 :        DO N = 1, nSpc
        1695           0 :           WRITE(MSG,*) TRIM(SpcNames(N)), ':', HcoIDs(N)
        1696           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        1697             :        ENDDO
        1698             :     ENDIF
        1699             :
        1700             :     ! Set up tracer indices
        1701           0 :     DO N = 1, nSpc
        1702           0 :        SELECT CASE( TRIM( SpcNames(N) ) )
        1703             :           CASE( 'POPG', 'POPG_BaP', 'POPG_PHE', 'POPG_PYR' )
        1704           0 :               Inst%IDTPOPG     = HcoIDs(N)
        1705             :           CASE( 'POPPOCPO', 'POPPOCPO_BaP', 'POPPOCPO_PHE', 'POPPOCPO_PYR' )
        1706           0 :               Inst%IDTPOPPOCPO = HcoIDs(N)
        1707             :           CASE( 'POPPBCPO', 'POPPBCPO_BaP', 'POPPBCPO_PHE', 'POPPBCPO_PYR' )
        1708           0 :              Inst%IDTPOPPBCPO = HcoIDs(N)
        1709             :           CASE DEFAULT
        1710             :              ! Do nothing
        1711             :        END SELECT
        1712             :     ENDDO
        1713             :
        1714             :     ! ERROR: POPG tracer is not found!
        1715           0 :     IF ( Inst%IDTPOPG <= 0 ) THEN
        1716           0 :        RC = HCO_FAIL
        1717           0 :        CALL HCO_ERROR( 'Cannot find POPG tracer in list of species!', RC )
        1718           0 :        RETURN
        1719             :     ENDIF
        1720             :
        1721             :     ! ERROR! POPPOCPO tracer is not found
        1722           0 :     IF ( Inst%IDTPOPPOCPO <= 0 ) THEN
        1723           0 :        RC = HCO_FAIL
        1724           0 :        CALL HCO_ERROR( 'Cannot find POPPOCPO tracer in list of species!', RC )
        1725           0 :        RETURN
        1726             :     ENDIF
        1727             :
        1728             :     ! ERROR! POPPBCPO tracer is not found
        1729           0 :     IF ( Inst%IDTPOPPBCPO <= 0 ) THEN
        1730           0 :        RC = HCO_FAIL
        1731           0 :        CALL HCO_ERROR( 'Cannot find POPPBCPO tracer in list of species!', RC )
        1732           0 :        RETURN
        1733             :     ENDIF
        1734             :
        1735             :     !=======================================================================
        1736             :     ! Activate this module and the fields of ExtState that it uses
        1737             :     !=======================================================================
        1738             :
        1739             :     ! Activate met fields required by this extension
        1740           0 :     ExtState%POPG%DoUse        = .TRUE.
        1741           0 :     ExtState%AIRVOL%DoUse      = .TRUE.
        1742           0 :     ExtState%AIRDEN%DoUse      = .TRUE.
        1743           0 :     ExtState%FRAC_OF_PBL%DoUse = .TRUE.
        1744           0 :     ExtState%LAI%DoUse         = .TRUE.
        1745           0 :     ExtState%PSC2_WET%DoUse    = .TRUE.
        1746           0 :     ExtState%SNOWHGT%DoUse     = .TRUE.
        1747           0 :     ExtState%TK%DoUse          = .TRUE.
        1748           0 :     ExtState%TSKIN%DoUse       = .TRUE.
        1749           0 :     ExtState%U10M%DoUse        = .TRUE.
        1750           0 :     ExtState%V10M%DoUse        = .TRUE.
        1751           0 :     ExtState%FRLAND%DoUse      = .TRUE.
        1752           0 :     ExtState%FRLANDIC%DoUse    = .TRUE.
        1753           0 :     ExtState%FROCEAN%DoUse     = .TRUE.
        1754           0 :     ExtState%FRSEAICE%DoUse    = .TRUE.
        1755           0 :     ExtState%FRLAKE%DoUse      = .TRUE.
        1756             :
        1757             :     !=======================================================================
        1758             :     ! Initialize data arrays
        1759             :     !=======================================================================
        1760             :
        1761           0 :     ALLOCATE( Inst%EmisPOPG( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
        1762           0 :     IF ( RC /= 0 ) THEN
        1763           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPG', RC )
        1764           0 :        RETURN
        1765             :     ENDIF
        1766           0 :     Inst%EmisPOPG = 0.0e0_hp
        1767             :
        1768           0 :     ALLOCATE( Inst%EmisPOPPOCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
        1769           0 :     IF ( RC /= 0 ) THEN
        1770           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPPOCPO', RC )
        1771           0 :        RETURN
        1772             :     ENDIF
        1773           0 :     Inst%EmisPOPPOCPO = 0.0e0_hp
        1774             :
        1775           0 :     ALLOCATE( Inst%EmisPOPPBCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
        1776           0 :     IF ( RC /= 0 ) THEN
        1777           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPPBCPO', RC )
        1778           0 :        RETURN
        1779             :     ENDIF
        1780           0 :     Inst%EmisPOPPBCPO = 0.0e0_hp
        1781             :
        1782           0 :     ALLOCATE( Inst%EmisPOPGfromSoil( HcoState%NX, HcoState%NY ), STAT=RC )
        1783           0 :     IF ( RC /= 0 ) THEN
        1784           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSoil', RC )
        1785           0 :        RETURN
        1786             :     ENDIF
        1787           0 :     Inst%EmisPOPGfromSoil = 0.0e0_hp
        1788             :
        1789           0 :     ALLOCATE( Inst%EmisPOPGfromLake( HcoState%NX, HcoState%NY ), STAT=RC )
        1790           0 :     IF ( RC /= 0 ) THEN
        1791           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLake', RC )
        1792           0 :        RETURN
        1793             :     ENDIF
        1794           0 :     Inst%EmisPOPGfromLake = 0.0e0_hp
        1795             :
        1796           0 :     ALLOCATE( Inst%EmisPOPGfromLeaf( HcoState%NX, HcoState%NY ), STAT=RC )
        1797           0 :     IF ( RC /= 0 ) THEN
        1798           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLeaf', RC )
        1799           0 :        RETURN
        1800             :     ENDIF
        1801           0 :     Inst%EmisPOPGfromLeaf = 0.0e0_hp
        1802             :
        1803           0 :     ALLOCATE( Inst%EmisPOPGfromSnow( HcoState%NX, HcoState%NY ), STAT=RC )
        1804           0 :     IF ( RC /= 0 ) THEN
        1805           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSnow', RC )
        1806           0 :        RETURN
        1807             :     ENDIF
        1808           0 :     Inst%EmisPOPGfromSnow = 0.0e0_hp
        1809             :
        1810           0 :     ALLOCATE( Inst%EmisPOPGfromOcean( HcoState%NX, HcoState%NY ), STAT=RC )
        1811           0 :     IF ( RC /= 0 ) THEN
        1812           0 :        CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromOcean', RC )
        1813           0 :        RETURN
        1814             :     ENDIF
        1815           0 :     Inst%EmisPOPGfromOcean = 0.0e0_hp
        1816             :
        1817           0 :     ALLOCATE( Inst%FluxPOPGfromSoilToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1818           0 :     IF ( RC /= 0 ) THEN
        1819           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromSoilToAir', RC )
        1820           0 :        RETURN
        1821             :     ENDIF
        1822           0 :     Inst%FluxPOPGfromSoilToAir = 0.0e0_hp
        1823             :
        1824           0 :     ALLOCATE( Inst%FluxPOPGfromAirToSoil( HcoState%NX, HcoState%NY ), STAT=RC )
        1825           0 :     IF ( RC /= 0 ) THEN
        1826           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToSoil', RC )
        1827           0 :        RETURN
        1828             :     ENDIF
        1829           0 :     Inst%FluxPOPGfromAirToSoil = 0.0e0_hp
        1830             :
        1831           0 :     ALLOCATE( Inst%FluxPOPGfromLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1832           0 :     IF ( RC /= 0 ) THEN
        1833           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLakeToAir', RC )
        1834           0 :        RETURN
        1835             :     ENDIF
        1836           0 :     Inst%FluxPOPGfromLakeToAir = 0.0e0_hp
        1837             :
        1838           0 :     ALLOCATE( Inst%FluxPOPGfromAirToLake( HcoState%NX, HcoState%NY ), STAT=RC )
        1839           0 :     IF ( RC /= 0 ) THEN
        1840           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLake', RC )
        1841           0 :        RETURN
        1842             :     ENDIF
        1843           0 :     Inst%FluxPOPGfromAirToLake = 0.0e0_hp
        1844             :
        1845           0 :     ALLOCATE( Inst%FluxPOPGfromLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1846           0 :     IF ( RC /= 0 ) THEN
        1847           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLeafToAir', RC )
        1848           0 :        RETURN
        1849             :     ENDIF
        1850           0 :     Inst%FluxPOPGfromLeafToAir = 0.0e0_hp
        1851             :
        1852           0 :     ALLOCATE( Inst%FluxPOPGfromAirToLeaf( HcoState%NX, HcoState%NY ), STAT=RC )
        1853           0 :     IF ( RC /= 0 ) THEN
        1854           0 :        CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLeaf', RC )
        1855           0 :        RETURN
        1856             :     ENDIF
        1857           0 :     Inst%FluxPOPGfromAirToLeaf = 0.0e0_hp
        1858             :
        1859           0 :     ALLOCATE( Inst%FugacitySoilToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1860           0 :     IF ( RC /= 0 ) THEN
        1861           0 :        CALL HCO_ERROR ( 'Cannot allocate FugacitySoilToAir', RC )
        1862           0 :        RETURN
        1863             :     ENDIF
        1864           0 :     Inst%FugacitySoilToAir = 0.0e0_hp
        1865             :
        1866           0 :     ALLOCATE( Inst%FugacityLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1867           0 :     IF ( RC /= 0 ) THEN
        1868           0 :        CALL HCO_ERROR ( 'Cannot allocate FugacityLakeToAir', RC )
        1869           0 :        RETURN
        1870             :     ENDIF
        1871           0 :     Inst%FugacityLakeToAir = 0.0e0_hp
        1872             :
        1873           0 :     ALLOCATE( Inst%FugacityLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC )
        1874           0 :     IF ( RC /= 0 ) THEN
        1875           0 :        CALL HCO_ERROR ( 'Cannot allocate FugacityLeafToAir', RC )
        1876           0 :        RETURN
        1877             :     ENDIF
        1878           0 :     Inst%FugacityLeafToAir = 0.0e0_hp
        1879             :
        1880             :     !=======================================================================
        1881             :     ! Create manual diagnostics
        1882             :     !=======================================================================
        1883           0 :     DO I = 1,12
        1884             :
        1885             :        ! Define diagnostic names. These have to match the names
        1886             :        ! in module HEMCO/Extensions/hcox_gc_POPs_mod.F90.
        1887             :        IF ( I == 1 ) THEN
        1888           0 :           DiagnName = 'EmisPOPGfromSoil'
        1889           0 :           OutUnit   = 'kg/m2/s'
        1890           0 :           Target2D  => Inst%EmisPOPGfromSoil
        1891             :        ELSEIF ( I == 2  ) THEN
        1892           0 :           DiagnName = 'EmisPOPGfromLake'
        1893           0 :           OutUnit   = 'kg/m2/s'
        1894           0 :           Target2D  => Inst%EmisPOPGfromLake
        1895             :        ELSEIF ( I == 3  ) THEN
        1896           0 :           DiagnName = 'EmisPOPGfromLeaf'
        1897           0 :           OutUnit   = 'kg/m2/s'
        1898           0 :           Target2D  => Inst%EmisPOPGfromLeaf
        1899             :        ELSEIF ( I == 4  ) THEN
        1900           0 :           DiagnName = 'FluxPOPGfromSoilToAir'
        1901           0 :           OutUnit   = 'ng/m2/day'
        1902           0 :           Target2D  => Inst%FluxPOPGfromSoilToAir
        1903             :        ELSEIF ( I == 5  ) THEN
        1904           0 :           DiagnName = 'FluxPOPGfromAirToSoil'
        1905           0 :           OutUnit   = 'ng/m2/day'
        1906           0 :           Target2D  => Inst%FluxPOPGfromAirToSoil
        1907             :        ELSEIF ( I == 6  ) THEN
        1908           0 :           DiagnName = 'FluxPOPGfromLakeToAir'
        1909           0 :           OutUnit   = 'ng/m2/day'
        1910           0 :           Target2D  => Inst%FluxPOPGfromLakeToAir
        1911             :        ELSEIF ( I == 7  ) THEN
        1912           0 :           DiagnName = 'FluxPOPGfromAirToLake'
        1913           0 :           OutUnit   = 'ng/m2/day'
        1914           0 :           Target2D  => Inst%FluxPOPGfromAirToLake
        1915             :        ELSEIF ( I == 8  ) THEN
        1916           0 :           DiagnName = 'FluxPOPGfromLeafToAir'
        1917           0 :           OutUnit   = 'ng/m2/day'
        1918           0 :           Target2D  => Inst%FluxPOPGfromLeafToAir
        1919             :        ELSEIF ( I == 9  ) THEN
        1920           0 :           DiagnName = 'FluxPOPGfromAirToLeaf'
        1921           0 :           OutUnit   = 'ng/m2/day'
        1922           0 :           Target2D  => Inst%FluxPOPGfromAirToLeaf
        1923             :        ELSEIF ( I == 10 ) THEN
        1924           0 :           DiagnName = 'FugacitySoilToAir'
        1925           0 :           OutUnit   = '1'
        1926           0 :           Target2D  => Inst%FugacitySoilToAir
        1927             :        ELSEIF ( I == 11 ) THEN
        1928           0 :           DiagnName = 'FugacityLakeToAir'
        1929           0 :           OutUnit   = '1'
        1930           0 :           Target2D  => Inst%FugacityLakeToAir
        1931             :        ELSEIF ( I == 12 ) THEN
        1932           0 :           DiagnName = 'FugacityLeafToAir'
        1933           0 :           OutUnit   = '1'
        1934           0 :           Target2D  => Inst%FugacityLeafToAir
        1935             :        ENDIF
        1936             :
        1937             :        ! Create manual diagnostics
        1938             :        CALL Diagn_Create( HcoState  = HcoState,          &
        1939             :                           cName     = TRIM( DiagnName ), &
        1940             :                           ExtNr     = ExtNr,             &
        1941             :                           Cat       = -1,                &
        1942             :                           Hier      = -1,                &
        1943             :                           HcoID     = -1,                &
        1944             :                           SpaceDim  = 2,                 &
        1945             :                           OutUnit   = OutUnit,           &
        1946             :                           AutoFill  = 0,                 &
        1947             :                           Trgt2D    = Target2D,          &
        1948           0 :                           RC        = RC )
        1949           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1950           0 :            CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
        1951           0 :            RETURN
        1952             :        ENDIF
        1953           0 :        Target2D => NULL()
        1954             :
        1955             :     ENDDO
        1956             :
        1957             :     !=======================================================================
        1958             :     ! Leave w/ success
        1959             :     !=======================================================================
        1960           0 :     IF ( ALLOCATED( HcoIDs   ) ) DEALLOCATE( HcoIDs   )
        1961           0 :     IF ( ALLOCATED( SpcNames ) ) DEALLOCATE( SpcNames )
        1962             :
        1963             :     ! Nullify pointers
        1964           0 :     Inst           => NULL()
        1965             :
        1966           0 :     CALL HCO_LEAVE( HcoState%Config%Err,RC )
        1967             :
        1968           0 :   END SUBROUTINE HCOX_GC_POPs_Init
        1969             : !EOC
        1970             : !------------------------------------------------------------------------------
        1971             : !                   Harmonized Emissions Component (HEMCO)                    !
        1972             : !------------------------------------------------------------------------------
        1973             : !BOP
        1974             : !
        1975             : ! !IROUTINE: HCOX_GC_POPs_Final
        1976             : !
        1977             : ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Final finalizes the HEMCO
        1978             : !  extension for the GEOS-Chem POPs specialty simulation.  All module
        1979             : !  arrays will be deallocated.
        1980             : !\\
        1981             : !\\
        1982             : ! !INTERFACE:
        1983             : !
        1984           0 :   SUBROUTINE HCOX_GC_POPs_Final( ExtState )
        1985             : !
        1986             : ! !INPUT PARAMETERS:
        1987             : !
        1988             :     TYPE(Ext_State),  POINTER       :: ExtState   ! Module options
        1989             : !
        1990             : ! !REVISION HISTORY:
        1991             : !  19 Aug 2014 - M. Sulprizio- Initial version
        1992             : !  See https://github.com/geoschem/hemco for complete history
        1993             : !EOP
        1994             : !------------------------------------------------------------------------------
        1995             : !BOC
        1996             :
        1997             :     !=======================================================================
        1998             :     ! HCOX_GC_POPs_FINAL begins here!
        1999             :     !=======================================================================
        2000             :
        2001           0 :     CALL InstRemove( ExtState%GC_POPs )
        2002             :
        2003           0 :   END SUBROUTINE HCOX_GC_POPs_Final
        2004             : !EOC
        2005             : !------------------------------------------------------------------------------
        2006             : !                   Harmonized Emissions Component (HEMCO)                    !
        2007             : !------------------------------------------------------------------------------
        2008             : !BOP
        2009             : !
        2010             : ! !IROUTINE: InstGet
        2011             : !
        2012             : ! !DESCRIPTION: Subroutine InstGet returns a poiner to the desired instance.
        2013             : !\\
        2014             : !\\
        2015             : ! !INTERFACE:
        2016             : !
        2017           0 :   SUBROUTINE InstGet ( Instance, Inst, RC, PrevInst )
        2018             : !
        2019             : ! !INPUT PARAMETERS:
        2020             : !
        2021             :     INTEGER                             :: Instance
        2022             :     TYPE(MyInst),     POINTER           :: Inst
        2023             :     INTEGER                             :: RC
        2024             :     TYPE(MyInst),     POINTER, OPTIONAL :: PrevInst
        2025             : !
        2026             : ! !REVISION HISTORY:
        2027             : !  18 Feb 2016 - C. Keller   - Initial version
        2028             : !  See https://github.com/geoschem/hemco for complete history
        2029             : !EOP
        2030             : !------------------------------------------------------------------------------
        2031             : !BOC
        2032             :     TYPE(MyInst),     POINTER    :: PrvInst
        2033             :
        2034             :     !=================================================================
        2035             :     ! InstGet begins here!
        2036             :     !=================================================================
        2037             :
        2038             :     ! Get instance. Also archive previous instance.
        2039           0 :     PrvInst => NULL()
        2040           0 :     Inst    => AllInst
        2041           0 :     DO WHILE ( ASSOCIATED(Inst) )
        2042           0 :        IF ( Inst%Instance == Instance ) EXIT
        2043           0 :        PrvInst => Inst
        2044           0 :        Inst    => Inst%NextInst
        2045             :     END DO
        2046           0 :     IF ( .NOT. ASSOCIATED( Inst ) ) THEN
        2047           0 :        RC = HCO_FAIL
        2048           0 :        RETURN
        2049             :     ENDIF
        2050             :
        2051             :     ! Pass output arguments
        2052           0 :     IF ( PRESENT(PrevInst) ) PrevInst => PrvInst
        2053             :
        2054             :     ! Cleanup & Return
        2055           0 :     PrvInst => NULL()
        2056           0 :     RC = HCO_SUCCESS
        2057             :
        2058             :   END SUBROUTINE InstGet
        2059             : !EOC
        2060             : !------------------------------------------------------------------------------
        2061             : !                   Harmonized Emissions Component (HEMCO)                    !
        2062             : !------------------------------------------------------------------------------
        2063             : !BOP
        2064             : !
        2065             : ! !IROUTINE: InstCreate
        2066             : !
        2067             : ! !DESCRIPTION: Subroutine InstCreate creates a new instance.
        2068             : !\\
        2069             : !\\
        2070             : ! !INTERFACE:
        2071             : !
        2072           0 :   SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC )
        2073             : !
        2074             : ! !INPUT PARAMETERS:
        2075             : !
        2076             :     INTEGER,       INTENT(IN)       :: ExtNr
        2077             : !
        2078             : ! !OUTPUT PARAMETERS:
        2079             : !
        2080             :     INTEGER,       INTENT(  OUT)    :: Instance
        2081             :     TYPE(MyInst),  POINTER          :: Inst
        2082             : !
        2083             : ! !INPUT/OUTPUT PARAMETERS:
        2084             : !
        2085             :     INTEGER,       INTENT(INOUT)    :: RC
        2086             : !
        2087             : ! !REVISION HISTORY:
        2088             : !  18 Feb 2016 - C. Keller   - Initial version
        2089             : !  See https://github.com/geoschem/hemco for complete history
        2090             : !EOP
        2091             : !------------------------------------------------------------------------------
        2092             : !BOC
        2093             :     TYPE(MyInst), POINTER          :: TmpInst
        2094             :     INTEGER                        :: nnInst
        2095             :
        2096             :     !=================================================================
        2097             :     ! InstCreate begins here!
        2098             :     !=================================================================
        2099             :
        2100             :     ! ----------------------------------------------------------------
        2101             :     ! Generic instance initialization
        2102             :     ! ----------------------------------------------------------------
        2103             :
        2104             :     ! Initialize
        2105           0 :     Inst => NULL()
        2106             :
        2107             :     ! Get number of already existing instances
        2108           0 :     TmpInst => AllInst
        2109           0 :     nnInst = 0
        2110           0 :     DO WHILE ( ASSOCIATED(TmpInst) )
        2111           0 :        nnInst  =  nnInst + 1
        2112           0 :        TmpInst => TmpInst%NextInst
        2113             :     END DO
        2114             :
        2115             :     ! Create new instance
        2116           0 :     ALLOCATE(Inst)
        2117           0 :     Inst%Instance = nnInst + 1
        2118           0 :     Inst%ExtNr    = ExtNr
        2119             :
        2120             :     ! Attach to instance list
        2121           0 :     Inst%NextInst => AllInst
        2122           0 :     AllInst       => Inst
        2123             :
        2124             :     ! Update output instance
        2125           0 :     Instance = Inst%Instance
        2126             :
        2127             :     ! ----------------------------------------------------------------
        2128             :     ! Type specific initialization statements follow below
        2129             :     ! ----------------------------------------------------------------
        2130             :
        2131             :     ! Return w/ success
        2132           0 :     RC = HCO_SUCCESS
        2133             :
        2134           0 :   END SUBROUTINE InstCreate
        2135             : !EOC
        2136             : !------------------------------------------------------------------------------
        2137             : !                   Harmonized Emissions Component (HEMCO)                    !
        2138             : !------------------------------------------------------------------------------
        2139             : !BOP
        2140             : !
        2141             : ! !IROUTINE: InstRemove
        2142             : !
        2143             : ! !DESCRIPTION: Subroutine InstRemove creates a new instance.
        2144             : !\\
        2145             : !\\
        2146             : ! !INTERFACE:
        2147             : !
        2148           0 :   SUBROUTINE InstRemove ( Instance )
        2149             : !
        2150             : ! !INPUT PARAMETERS:
        2151             : !
        2152             :     INTEGER                         :: Instance
        2153             : !
        2154             : ! !REVISION HISTORY:
        2155             : !  18 Feb 2016 - C. Keller   - Initial version
        2156             : !  See https://github.com/geoschem/hemco for complete history
        2157             : !EOP
        2158             : !------------------------------------------------------------------------------
        2159             : !BOC
        2160             :     INTEGER                     :: RC
        2161             :     TYPE(MyInst), POINTER       :: PrevInst
        2162             :     TYPE(MyInst), POINTER       :: Inst
        2163             :
        2164             :     !=================================================================
        2165             :     ! InstRemove begins here!
        2166             :     !=================================================================
        2167             :
        2168             :     ! Init
        2169           0 :     PrevInst => NULL()
        2170           0 :     Inst     => NULL()
        2171             :
        2172             :     ! Get instance. Also archive previous instance.
        2173           0 :     CALL InstGet ( Instance, Inst, RC, PrevInst=PrevInst )
        2174             :
        2175             :     ! Instance-specific deallocation
        2176           0 :     IF ( ASSOCIATED(Inst) ) THEN
        2177             :
        2178             :        !---------------------------------------------------------------------
        2179             :        ! Deallocate fields of Inst before popping off from the list
        2180             :        ! in order to avoid memory leaks (Bob Yantosca (17 Aug 2022)
        2181             :        !---------------------------------------------------------------------
        2182           0 :        IF ( ASSOCIATED( Inst%EmisPOPPOCPO ) ) THEN
        2183           0 :           DEALLOCATE( Inst%EmisPOPPOCPO )
        2184             :        ENDIF
        2185           0 :        Inst%EmisPOPPOCPO => NULL()
        2186             :
        2187           0 :        IF ( ASSOCIATED( Inst%EmisPOPPBCPO ) ) THEN
        2188           0 :           DEALLOCATE( Inst%EmisPOPPBCPO )
        2189             :        ENDIF
        2190           0 :        Inst%EmisPOPPBCPO => NULL()
        2191             :
        2192           0 :        IF ( ASSOCIATED( Inst%EmisPOPG ) ) THEN
        2193           0 :           DEALLOCATE( Inst%EmisPOPG )
        2194             :        ENDIF
        2195           0 :        Inst%EmisPOPG => NULL()
        2196             :
        2197           0 :        IF ( ASSOCIATED( Inst%EmisPOPGfromSoil ) ) THEN
        2198           0 :           DEALLOCATE( Inst%EmisPOPGfromSoil )
        2199             :        ENDIF
        2200           0 :        Inst%EmisPOPGfromSoil => NULL()
        2201             :
        2202           0 :        IF ( ASSOCIATED( Inst%EmisPOPGfromLake ) ) THEN
        2203           0 :           DEALLOCATE( Inst%EmisPOPGfromLake )
        2204             :        ENDIF
        2205           0 :        Inst%EmisPOPGfromLake => NULL()
        2206             :
        2207           0 :        IF ( ASSOCIATED( Inst%EmisPOPGfromLeaf ) ) THEN
        2208           0 :           DEALLOCATE( Inst%EmisPOPGfromLeaf )
        2209             :        ENDIF
        2210           0 :        Inst%EmisPOPGfromLeaf => NULL()
        2211             :
        2212           0 :        IF ( ASSOCIATED( Inst%EmisPOPGfromSnow  ) ) THEN
        2213           0 :           DEALLOCATE( Inst%EmisPOPGfromSnow  )
        2214             :        ENDIF
        2215           0 :        Inst%EmisPOPGfromSnow => NULL()
        2216             :
        2217           0 :        IF ( ASSOCIATED( Inst%EmisPOPGfromOcean ) ) THEN
        2218           0 :           DEALLOCATE( Inst%EmisPOPGfromOcean    )
        2219             :        ENDIF
        2220           0 :        Inst%EmisPOPGfromOcean => NULL()
        2221             :
        2222           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromSoilToAir ) ) THEN
        2223           0 :           DEALLOCATE( Inst%FluxPOPGfromSoilToAir )
        2224             :        ENDIF
        2225           0 :        Inst%FluxPOPGfromSoilToAir => NULL()
        2226             :
        2227           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromAirToSoil ) ) THEN
        2228           0 :           DEALLOCATE( Inst%FluxPOPGfromAirToSoil )
        2229             :        ENDIF
        2230           0 :        Inst%FluxPOPGfromAirToSoil => NULL()
        2231             :
        2232           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromLakeToAir ) ) THEN
        2233           0 :           DEALLOCATE( Inst%FluxPOPGfromLakeToAir )
        2234             :        ENDIF
        2235           0 :        Inst%FluxPOPGfromLakeToAir => NULL()
        2236             :
        2237           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromAirToLake ) ) THEN
        2238           0 :           DEALLOCATE( Inst%FluxPOPGfromAirToLake )
        2239             :        ENDIF
        2240           0 :        Inst%FluxPOPGfromAirToLake => NULL()
        2241             :
        2242           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromLeafToAir ) ) THEN
        2243           0 :           DEALLOCATE( Inst%FluxPOPGfromLeafToAir )
        2244             :        ENDIF
        2245           0 :        Inst%FluxPOPGfromLeafToAir => NULL()
        2246             :
        2247           0 :        IF ( ASSOCIATED( Inst%FluxPOPGfromAirtoLeaf ) ) THEN
        2248           0 :           DEALLOCATE( Inst%FluxPOPGfromAirtoLeaf )
        2249             :        ENDIF
        2250           0 :        Inst%FluxPOPGfromAirtoLeaf => NULL()
        2251             :
        2252           0 :        IF ( ASSOCIATED( Inst%FugacitySoilToAir ) ) THEN
        2253           0 :           DEALLOCATE( Inst%FugacitySoilToAir )
        2254             :        ENDIF
        2255           0 :        Inst%FugacitySoilToAir => NULL()
        2256             :
        2257           0 :        IF ( ASSOCIATED( Inst%FugacityLakeToAir ) ) THEN
        2258           0 :           DEALLOCATE( Inst%FugacityLakeToAir )
        2259             :        ENDIF
        2260           0 :        Inst%FugacityLakeToAir => NULL()
        2261             :
        2262           0 :        IF ( ASSOCIATED( Inst%FugacityLeafToAir ) ) THEN
        2263           0 :           DEALLOCATE( Inst%FugacityLeafToAir )
        2264             :        ENDIF
        2265           0 :        Inst%FugacityLeafToAir => NULL()
        2266             :
        2267             :        !---------------------------------------------------------------------
        2268             :        ! Pop off instance from list
        2269             :        !---------------------------------------------------------------------
        2270           0 :        IF ( ASSOCIATED(PrevInst) ) THEN
        2271           0 :           PrevInst%NextInst => Inst%NextInst
        2272             :        ELSE
        2273           0 :           AllInst => Inst%NextInst
        2274             :        ENDIF
        2275           0 :        DEALLOCATE(Inst)
        2276             :     ENDIF
        2277             :
        2278             :     ! Free pointers before exiting
        2279           0 :     PrevInst => NULL()
        2280           0 :     Inst     => NULL()
        2281             :
        2282           0 :    END SUBROUTINE InstRemove
        2283             : !EOC
        2284           0 : END MODULE HCOX_GC_POPs_Mod

Generated by: LCOV version 1.14