Difference between revisions of "POPs simulation"

From Geos-chem
Jump to: navigation, search
(Prevent error when reading global OC)
(Previous issues that are now resolved)
Line 46: Line 46:
  
 
--[[User:Melissa Payer|Melissa Sulprizio]] 10:33, 24 February 2014 (EST)
 
--[[User:Melissa Payer|Melissa Sulprizio]] 10:33, 24 February 2014 (EST)
 +
 +
== Outstanding issues ==
 +
 +
=== Avoid div-by-zero errors in pops_mod.F ===
 +
 +
While testing the POPs simulation with the [[GEOS-Chem Unit Tester]], we encountered some division-by-zero errors.  The fixes are as follows:
 +
 +
(1) At about line 1661 of routine <tt>EMISSPOPS</tt> (in module file <tt>GeosCore/pops_mod.F</tt>), change this line of code:
 +
 +
          ! Check that sum thru PBL is equal to original emissions array
 +
          SUM_OF_ALL(I,J) = POP_TOT_EM(I,J) / SUM_OF_ALL(I,J)
 +
 +
to this:
 +
 +
          ! Check that sum thru PBL is equal to original emissions array
 +
          ! NOTE: Prevent div-by-zero floating point error (bmy, 4/14/14)
 +
          IF ( SUM_OF_ALL(I,J) > 0d0 ) THEN
 +
            SUM_OF_ALL(I,J) = POP_TOT_EM(I,J) / SUM_OF_ALL(I,J)
 +
          ENDIF
 +
 +
This prevents us from dividing by SUM_OF_ALL(I,J) if it’s zero.  Otherwise some kind of junk value or NaN could possibly propagate thru the code.
 +
 +
 +
(2) In routine <tt>CHEM_POPGP</tt> (in module GeosCore/pops_mod.F), make the following modifications:
 +
 +
      USE ERROR_MOD,          ONLY : SAFE_DIV
 +
 +
      . . .
 +
 +
      REAL*8          :: DENOM
 +
 +
      . . .
 +
 +
          ! Get AIRVOL
 +
          AIR_VOL = State_Met%AIRVOL(I,J,L)
 +
 +
!-----------------------------------------------------------------------------
 +
! Prior to 4/14/14:
 +
! Need to put error traps to prevent div-by-zero (bmy, 4/14/14)
 +
!        ! Define volume ratios:
 +
!        ! VR_OC_AIR = volume ratio of OC to air [unitless]   
 +
!        VR_OC_AIR  = C_OC_CHEM1 / AIR_VOL ! could be zero
 +
!
 +
!        ! VR_OC_BC  = volume ratio of OC to BC [unitless]
 +
!        VR_OC_BC = C_OC_CHEM1 / C_BC_CHEM1 ! could be zero or undefined
 +
!
 +
!        ! VR_BC_AIR = volume ratio of BC to air [unitless]
 +
!        VR_BC_AIR  = VR_OC_AIR / VR_OC_BC ! could be zero or undefined
 +
!
 +
!        ! VR_BC_OC  = volume ratio of BC to OC [unitless]
 +
!        VR_BC_OC    = 1d0 / VR_OC_BC ! could be zero or undefined
 +
!
 +
!        ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase,
 +
!        ! BC-phase, and gas phase with new time step (should only change if
 +
!        ! temp changes or OC/BC concentrations change)
 +
!        OC_AIR_RATIO = 1d0 / (KOA_T * VR_OC_AIR)
 +
!        OC_BC_RATIO = 1d0 / (KOC_BC_T * VR_OC_BC)
 +
!
 +
!        BC_AIR_RATIO = 1d0 / (KBC_T * VR_BC_AIR)
 +
!        BC_OC_RATIO = 1d0 / (KBC_OC_T * VR_BC_OC)
 +
!-----------------------------------------------------------------------------
 +
 +
          ! Define volume ratios:
 +
          ! VR_OC_AIR = volume ratio of OC to air [unitless]   
 +
          VR_OC_AIR    = C_OC_CHEM1 / AIR_VOL ! could be zero
 +
 +
          ! VR_OC_BC = volume ratio of OC to BC [unitless]
 +
          VR_OC_BC      = SAFE_DIV( C_OC_CHEM1, C_BC_CHEM1, 0d0 )
 +
 +
          ! VR_BC_AIR = volume ratio of BC to air [unitless]
 +
          VR_BC_AIR    = SAFE_DIV( VR_OC_AIR,  VR_OC_BC,  0d0 )
 +
 +
          ! VR_BC_OC = volume ratio of BC to OC [unitless]
 +
          VR_BC_OC      = SAFE_DIV( 1d0,        VR_OC_BC,  0d0 )
 +
 +
          ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase,
 +
          ! BC-phase, and gas phase with new time step (should only change if
 +
          ! temp changes or OC/BC concentrations change)
 +
          DENOM        = KOA_T * VR_OC_AIR
 +
          OC_AIR_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )
 +
 +
          DENOM        = KOC_BC_T * VR_OC_BC
 +
          OC_BC_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )
 +
 +
          DENOM        = KBC_T * VR_BC_AIR
 +
          BC_AIR_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )
 +
 +
          DENOM        = KBC_OC_T * VR_BC_OC
 +
          BC_OC_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )
 +
 +
The function <tt>SAFE_DIV</tt> (in module file <tt>GeosUtil/error_mod.F</tt>) tests if the division can be done, and if not, it will assign an alternate value (in this case, 0d0).  I needed to add these in to prevent div-by-zero errors from happening in the computations for variables: <tt>VR_OC_AIR, VR_BC_AIR, VR_BC_OC, VR_OC_BC, OC_AIR_RATIO, OC_BC_RATIO, BC_AIR_RATIO, BC_OC_RATIO</tt>.
 +
 +
--[[User:Bmy|Bob Y.]] 14:07, 14 April 2014 (EDT)
  
 
== References ==
 
== References ==
  
 
#Friedman, C.L. and N.E. Selin. ''Long-range atmospheric transport of polycyclic aromatic hydrocarbons: A global 3-D model analysis including evaluation of Arctic sources'', <u>Environ. Sci. Technol.</u>, '''46''', 9501-9510, 2012.
 
#Friedman, C.L. and N.E. Selin. ''Long-range atmospheric transport of polycyclic aromatic hydrocarbons: A global 3-D model analysis including evaluation of Arctic sources'', <u>Environ. Sci. Technol.</u>, '''46''', 9501-9510, 2012.

Revision as of 18:07, 14 April 2014

On this page we include information relevant to the persistent organic pollutant (POPs) simulation in GEOS-Chem.

Overview

Authors and collaborators

PAH simulation

This update was tested in the 1-month benchmark simulation v9-02c and approved on 29 Nov 2012.

A simulation for polycyclic aromatic hydrocarbons (PAHs) has been added to GEOS-Chem v8-03-02 following Friedman and Selin (2012). Three separate PAHs are modeled: phenanthrene (PHE), pyrene (PYR), and benzo[a]pyrene (BaP). Separate modules have been created for each of the PAHs. Users can choose between different compounds in the input file.

--Melissa Payer 12:31, 21 September 2012 (EDT)

PCB simulation

Info to be added.

Previous issues that are now resolved

Prevent error when reading global OC

NOTE: This issue was resolved in the GEOS-Chem v9-02 public release (01 Mar 2014).

In routine GET_GLOBAL_OC (in module GeosCore/global_oc_mod.F), we have:

     ! Data is only available for 2005-2009
     IF ( THISYEAR < 2005 ) THEN
        YEAR = 2005
     ELSE IF ( THIS YEAR > 2009 ) THEN
        YEAR = 2009
     ELSE
        YEAR = THISYEAR
     ENDIF

     ! Get the TAU0 value for the start of the given month
     XTAU = GET_TAU0( THISMONTH, 1, THISYEAR )

To avoid a segmentation fault error when THISYEAR is outside of 2005-2009, we must change XTAU to:

     ! Get the TAU0 value for the start of the given month
     XTAU = GET_TAU0( THISMONTH, 1, YEAR )

--Melissa Sulprizio 10:33, 24 February 2014 (EST)

Outstanding issues

Avoid div-by-zero errors in pops_mod.F

While testing the POPs simulation with the GEOS-Chem Unit Tester, we encountered some division-by-zero errors. The fixes are as follows:

(1) At about line 1661 of routine EMISSPOPS (in module file GeosCore/pops_mod.F), change this line of code:

         ! Check that sum thru PBL is equal to original emissions array
         SUM_OF_ALL(I,J) = POP_TOT_EM(I,J) / SUM_OF_ALL(I,J)

to this:

         ! Check that sum thru PBL is equal to original emissions array
         ! NOTE: Prevent div-by-zero floating point error (bmy, 4/14/14)
         IF ( SUM_OF_ALL(I,J) > 0d0 ) THEN 
            SUM_OF_ALL(I,J) = POP_TOT_EM(I,J) / SUM_OF_ALL(I,J)
         ENDIF

This prevents us from dividing by SUM_OF_ALL(I,J) if it’s zero. Otherwise some kind of junk value or NaN could possibly propagate thru the code.


(2) In routine CHEM_POPGP (in module GeosCore/pops_mod.F), make the following modifications:

      USE ERROR_MOD,          ONLY : SAFE_DIV

      . . .

      REAL*8          :: DENOM
      . . .
         ! Get AIRVOL
         AIR_VOL = State_Met%AIRVOL(I,J,L)
!-----------------------------------------------------------------------------
! Prior to 4/14/14:
! Need to put error traps to prevent div-by-zero (bmy, 4/14/14)
!         ! Define volume ratios:
!         ! VR_OC_AIR = volume ratio of OC to air [unitless]     
!         VR_OC_AIR   = C_OC_CHEM1 / AIR_VOL ! could be zero
!
!         ! VR_OC_BC  = volume ratio of OC to BC [unitless]
!         VR_OC_BC = C_OC_CHEM1 / C_BC_CHEM1 ! could be zero or undefined
!
!         ! VR_BC_AIR = volume ratio of BC to air [unitless]
!         VR_BC_AIR   = VR_OC_AIR / VR_OC_BC ! could be zero or undefined
!
!         ! VR_BC_OC  = volume ratio of BC to OC [unitless]
!         VR_BC_OC    = 1d0 / VR_OC_BC ! could be zero or undefined
!
!         ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase, 
!         ! BC-phase, and gas phase with new time step (should only change if 
!         ! temp changes or OC/BC concentrations change) 
!         OC_AIR_RATIO = 1d0 / (KOA_T * VR_OC_AIR) 
!         OC_BC_RATIO = 1d0 / (KOC_BC_T * VR_OC_BC) 
!
!         BC_AIR_RATIO = 1d0 / (KBC_T * VR_BC_AIR) 
!         BC_OC_RATIO = 1d0 / (KBC_OC_T * VR_BC_OC)
!-----------------------------------------------------------------------------

         ! Define volume ratios:
         ! VR_OC_AIR = volume ratio of OC to air [unitless]     
         VR_OC_AIR     = C_OC_CHEM1 / AIR_VOL ! could be zero

         ! VR_OC_BC = volume ratio of OC to BC [unitless]
         VR_OC_BC      = SAFE_DIV( C_OC_CHEM1, C_BC_CHEM1, 0d0 )

         ! VR_BC_AIR = volume ratio of BC to air [unitless]
         VR_BC_AIR     = SAFE_DIV( VR_OC_AIR,  VR_OC_BC,   0d0 )

         ! VR_BC_OC = volume ratio of BC to OC [unitless]
         VR_BC_OC      = SAFE_DIV( 1d0,        VR_OC_BC,   0d0 ) 

         ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase, 
         ! BC-phase, and gas phase with new time step (should only change if 
         ! temp changes or OC/BC concentrations change) 
         DENOM         = KOA_T * VR_OC_AIR
         OC_AIR_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )
         DENOM         = KOC_BC_T * VR_OC_BC
         OC_BC_RATIO   = SAFE_DIV( 1d0,        DENOM,      0d0 )

         DENOM         = KBC_T * VR_BC_AIR
         BC_AIR_RATIO  = SAFE_DIV( 1d0,        DENOM,      0d0 )

         DENOM         = KBC_OC_T * VR_BC_OC
         BC_OC_RATIO   = SAFE_DIV( 1d0,        DENOM,      0d0 )

The function SAFE_DIV (in module file GeosUtil/error_mod.F) tests if the division can be done, and if not, it will assign an alternate value (in this case, 0d0). I needed to add these in to prevent div-by-zero errors from happening in the computations for variables: VR_OC_AIR, VR_BC_AIR, VR_BC_OC, VR_OC_BC, OC_AIR_RATIO, OC_BC_RATIO, BC_AIR_RATIO, BC_OC_RATIO.

--Bob Y. 14:07, 14 April 2014 (EDT)

References

  1. Friedman, C.L. and N.E. Selin. Long-range atmospheric transport of polycyclic aromatic hydrocarbons: A global 3-D model analysis including evaluation of Arctic sources, Environ. Sci. Technol., 46, 9501-9510, 2012.