Mercury: Difference between revisions

From Geos-chem
Jump to navigation Jump to search
No edit summary
Line 81: Line 81:


The only thing I have to add is that at first I didn't realize that wet deposition of Hg(II) was composed of both WETDCV and WETDLS.  Is it important to save those components out as 2 separate parts? (nvd)
The only thing I have to add is that at first I didn't realize that wet deposition of Hg(II) was composed of both WETDCV and WETDLS.  Is it important to save those components out as 2 separate parts? (nvd)
== Bug fixes ==
''''' [mailto:ccc@seas.harvard.edu Claire Carouge] wrote:'''''
:I came across a parallelization bug in the snowpack routine from Chris's code (in mercury_mod.f).
:The initial lines of code are:
      ! Emit Hg(0) at a steady rate, based on 180 d residence
      ! time in snowpack, based on cycle observed at Alert
      ! (e.g. Steffen et al. 2008)
      K_EMIT = 6D-8
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, NN )
!$OMP+PRIVATE( SNOW_HG_NEW, JLOOP, K_EMIT )
      DO J  = 1, JJPAR
      DO I  = 1, IIPAR
          ! 1-D grid box index for SUNCOS
          JLOOP = ( (J-1) * IIPAR ) + I
          ! If the sun is set, then no emissions, go to next box
          IF (SUNCOS(JLOOP)<0D0) CYCLE
          ! Decrease residence time to 1 week when T > -3C
          IF (T(I,J,1) > 270D0) K_EMIT = 1.6D-6
:But this creates problems for the value of K_EMIT for several reasons.
:1/ When declaring K_EMIT as PRIVATE, the code creates a variable of the same type as K_EMIT for each processor (thread) but doesn't give it any value! So if T(I,J,1) <= 270D0, then K_EMIT in the parallel loop has no value. The initialization before the loop is completely ignored.
:2/ Let's imagine one processor starts by handling the point (1,1) and you have T(1,1,1) > 270D0 then K_EMIT = 1.6D-6. Then the same processor might handle the point (1,2) and let's say T(1,2,1) <= 270D0, then I guess you would like to have K_EMIT = 6D-8. But there is no line in the J,I loop to change the value of K_EMIT back to 6D-8. So this processor keeps using 1.6D-6. The K_EMIT is private to each processor not each box.
:So to conclude, in this loop, K_EMIT has either no value or 1.6D-6. And sometimes K_EMIT = 1.6D-6 when it should be 6D-8.
:Here is the fix:
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, NN )
!$OMP+PRIVATE( SNOW_HG_NEW, JLOOP, K_EMIT )
      DO J  = 1, JJPAR
      DO I  = 1, IIPAR
          ! 1-D grid box index for SUNCOS
          JLOOP = ( (J-1) * IIPAR ) + I
          ! If the sun is set, then no emissions, go to next box
          IF (SUNCOS(JLOOP)<0D0) CYCLE
          ! Decrease residence time to 1 week when T > -3C
          IF (T(I,J,1) > 270D0) THEN
            K_EMIT = 1.6D-6
          ELSE
            ! Emit Hg(0) at a steady rate, based on 180 d residence
            ! time in snowpack, based on cycle observed at Alert
            ! (e.g. Steffen et al. 2008)
            K_EMIT = 6D-8
          ENDIF
:So the difference is the "ELSE" part for the condition on the temperature. And we can get rid of the initialization before the loop (useless).
:I haven't run for long simulations so I'm not sure how much impact this has on the results.
--[[User:Ccarouge|Ccarouge]] 15:50, 4 May 2010 (EDT)

Revision as of 19:50, 4 May 2010

Hello Team-Hg

Code

Here's the list of things I *think* we came up with on Friday regarding cleaning up the code (feel free to add/comment)

1) Standardize the solver. Everyone should use the solver Chris developed for the Hg chemistry
(located at ~cdh/GC/RevisedChem.v7-04-06/mercury_mod.f)

  • outstanding issue - dry dep of Hg0* <- currently working on this (eds)


2) Catalog all emissions options and develop clear flagging system to choose your own adventure. This will include:

  • different anthropogenic emissions scenarios/corrections (i.e Jaffe vs. Streets) <- going to work on the anthro emissions soon (eds),
  • different land emissions. <- nvd will work on this


3) Diagnostics. See separate section below.
4) Comment everything in the code. Remove old bits of code that are hanging around & commented out.
5) GEOS-5
6) Get the land stuff out of ocean_mercury_mod.f and into it's own module <- nvd will work on this

Chemistry 'issues'

1) Oxidant. Chris has a simulation with Hg-Br chemistry and SS aerosol deposition; the global budget is ok, but the Br concentrations in the BL are too low to generate diurnal cycles. (cdh working on it)
2) Snow/ice scavenging of HgII
3) Dry deposition of "aqueous HgII. (Explanation from Chris: We calculate the fraction of HgII, Fg, which is gas phase. But we're currently applying the dry deposition velocity to both gas and aqueous fractions. I think it would be better to deposit the aqueous fraction at the velocity of HgP; this would be slower dep, but I don't know how much. This is definitely up for discussion.)

Diagnostics

Here are some suggested changes:

  1. Emissions should have units 'kg/m2/s' or something of the form 'mass/area/time' (they are currently 'kg'). The HG-SRCE diagnostic currently has all of the Ocean tracers and fluxes; these should go elsewhere.
  2. The Ocean Hg0, Hg2, HgC should have concentration units not kg. Is 'molar' the best choice? Fluxes of these should be in concentration/time, not kg.
  3. The ocean restart files should have concentration units not kg. They currently use the category 'OCEAN-HG' which would make sense for the ND03 ocean Hg0, Hg2, HgC output too.
  4. 'PL-HG2-$' doesn't really describe all of the fluxes in our model. There are a lot of diagnostic quantities which are either chemical P/L fluxes or rate constants. I think these should all be in one diagnostic called something like 'PL-HG-$' (or maybe 'PL-HG-A', 'PL-HG-O' to separate the atmosphere and ocean). The fluxes in this diagnostic would include redox in air and water, colloidal sinking, ocean-atmosphere piston velocity, ...

Here are the current GEOS-Chem Hg outputs 

      CATEGORY ILUN TRCNAME   TRC         UNIT      TAU0(DATE)       DIMENSIONS
  1 : IJ-AVG-$   23     Hg0     1         pptv 157776.00(2003010100)  72 46 30
  2 : IJ-AVG-$   23     Hg2     2         pptv 157776.00(2003010100)  72 46 30
  3 : IJ-AVG-$   23     HgP     3         pptv 157776.00(2003010100)  72 46 30
  4 : WETDCV-$   23     Hg2  3002         kg/s 157776.00(2003010100)  72 46 30
  5 : WETDCV-$   23     HgP  3003         kg/s 157776.00(2003010100)  72 46 30
  6 : WETDLS-$   23     Hg2  3002         kg/s 157776.00(2003010100)  72 46 30
  7 : WETDLS-$   23     HgP  3003         kg/s 157776.00(2003010100)  72 46 30
  8 :  HG-SRCE   23  Hg0_an 34001           kg 157776.00(2003010100)  72 46  1
  9 :  HG-SRCE   23  Hg0_aq 34002           kg 157776.00(2003010100)  72 46  1
 10 :  HG-SRCE   23  Hg0_oc 34003           kg 157776.00(2003010100)  72 46  1
 11 :  HG-SRCE   23  Hg0_ln 34004           kg 157776.00(2003010100)  72 46  1
 12 :  HG-SRCE   23  Hg0_na 34005           kg 157776.00(2003010100)  72 46  1
 13 :  HG-SRCE   23  Hg2_an 34006           kg 157776.00(2003010100)  72 46  1
 14 :  HG-SRCE   23  Hg2_aq 34007           kg 157776.00(2003010100)  72 46  1
 15 :  HG-SRCE   23  Hg2_sk 34008           kg 157776.00(2003010100)  72 46  1
 16 :  HG-SRCE   23  HgP_an 34009           kg 157776.00(2003010100)  72 46  1
 17 :  HG-SRCE   23    KwHg 34010         cm/h 157776.00(2003010100)  72 46  1
 18 :  HG-SRCE   23     HgC 34011           kg 157776.00(2003010100)  72 46  1
 19 :  HG-SRCE   23 Hg_to_C 34012           kg 157776.00(2003010100)  72 46  1
 20 : PL-HG2-$   23 Hg2_Hg0 35001           kg 157776.00(2003010100)  72 46 30
 21 : PL-HG2-$   23  Hg2_OH 35002           kg 157776.00(2003010100)  72 46 30
 22 : PL-HG2-$   23  Hg2_O3 35003           kg 157776.00(2003010100)  72 46 30
 23 : PL-HG2-$   23  Hg2_SS 35004           kg 157776.00(2003010100)  72 46  1
 24 : PL-HG2-$   23 Hg2_SSR 35005           /s 157776.00(2003010100)  72 46  1
 25 : DRYD-FLX   23   Hg0df 36001  molec/cm2/s 157776.00(2003010100)  72 46  1
 26 : DRYD-FLX   23   Hg2df 36002  molec/cm2/s 157776.00(2003010100)  72 46  1
 27 : DRYD-FLX   23   HgPdf 36003  molec/cm2/s 157776.00(2003010100)  72 46  1
 28 : DRYD-VEL   23   Hg2dv 37002         cm/s 157776.00(2003010100)  72 46  1
 29 : DRYD-VEL   23   HgPdv 37003         cm/s 157776.00(2003010100)  72 46  1

I think we should change lines 8-24 (I've kept the same line numbers and TRCNAME, but changed CATEGORY, TRC, or UNIT): 

      CATEGORY  TRCNAME    TRC         UNIT 
  8 :  HG-SRCE   Hg0_an  34001      kg/m2/s 
 10 :  HG-SRCE   Hg0_oc  34002      kg/m2/s 
 11 :  HG-SRCE   Hg0_ln  34003      kg/m2/s 
 12 :  HG-SRCE   Hg0_na  34004      kg/m2/s 
 13 :  HG-SRCE   Hg2_an  34005      kg/m2/s
 16 :  HG-SRCE   HgP_an  34006      kg/m2/s 
  9 : OCEAN-HG   Hg0_aq  xxxx1        mol/L 
 14 : OCEAN-HG   Hg2_aq  xxxx2        mol/L 
 18 : OCEAN-HG   HgC     xxxx3        mol/L 
 20 :  PL-HG-A   Hg2_Hg0 35001      kg/m3/s 
 21 :  PL-HG-A   Hg2_OH  35002      kg/m3/s 
 22 :  PL-HG-A   Hg2_O3  35003      kg/m3/s 
 23 :  PL-HG-A   Hg2_SS  35004      kg/m3/s 
 24 :  PL-HG-A   Hg2_SSR 35005           /s
 15 :  PL-HG-O   Hg2_sk  xxxx1      kg/m3/s 
 19 :  PL-HG-O   Hg_to_C xxxx2      kg/m3/s
 17 :  PL-HG-O   KwHg    xxxx3         cm/h


The only thing I have to add is that at first I didn't realize that wet deposition of Hg(II) was composed of both WETDCV and WETDLS. Is it important to save those components out as 2 separate parts? (nvd)


Bug fixes

Claire Carouge wrote:

I came across a parallelization bug in the snowpack routine from Chris's code (in mercury_mod.f).
The initial lines of code are:
      ! Emit Hg(0) at a steady rate, based on 180 d residence
      ! time in snowpack, based on cycle observed at Alert 
      ! (e.g. Steffen et al. 2008)
      K_EMIT = 6D-8

!$OMP PARALLEL DO 
!$OMP+DEFAULT( SHARED ) 
!$OMP+PRIVATE( I, J, NN )
!$OMP+PRIVATE( SNOW_HG_NEW, JLOOP, K_EMIT )
      DO J  = 1, JJPAR
      DO I  = 1, IIPAR

         ! 1-D grid box index for SUNCOS
         JLOOP = ( (J-1) * IIPAR ) + I

         ! If the sun is set, then no emissions, go to next box
         IF (SUNCOS(JLOOP)<0D0) CYCLE 

         ! Decrease residence time to 1 week when T > -3C
         IF (T(I,J,1) > 270D0) K_EMIT = 1.6D-6
But this creates problems for the value of K_EMIT for several reasons.
1/ When declaring K_EMIT as PRIVATE, the code creates a variable of the same type as K_EMIT for each processor (thread) but doesn't give it any value! So if T(I,J,1) <= 270D0, then K_EMIT in the parallel loop has no value. The initialization before the loop is completely ignored.
2/ Let's imagine one processor starts by handling the point (1,1) and you have T(1,1,1) > 270D0 then K_EMIT = 1.6D-6. Then the same processor might handle the point (1,2) and let's say T(1,2,1) <= 270D0, then I guess you would like to have K_EMIT = 6D-8. But there is no line in the J,I loop to change the value of K_EMIT back to 6D-8. So this processor keeps using 1.6D-6. The K_EMIT is private to each processor not each box.
So to conclude, in this loop, K_EMIT has either no value or 1.6D-6. And sometimes K_EMIT = 1.6D-6 when it should be 6D-8.
Here is the fix:
!$OMP PARALLEL DO 
!$OMP+DEFAULT( SHARED ) 
!$OMP+PRIVATE( I, J, NN )
!$OMP+PRIVATE( SNOW_HG_NEW, JLOOP, K_EMIT )
      DO J  = 1, JJPAR
      DO I  = 1, IIPAR

         ! 1-D grid box index for SUNCOS
         JLOOP = ( (J-1) * IIPAR ) + I

         ! If the sun is set, then no emissions, go to next box
         IF (SUNCOS(JLOOP)<0D0) CYCLE 

         ! Decrease residence time to 1 week when T > -3C
         IF (T(I,J,1) > 270D0) THEN
            K_EMIT = 1.6D-6
         ELSE
            ! Emit Hg(0) at a steady rate, based on 180 d residence
            ! time in snowpack, based on cycle observed at Alert 
            ! (e.g. Steffen et al. 2008)
            K_EMIT = 6D-8
         ENDIF
So the difference is the "ELSE" part for the condition on the temperature. And we can get rid of the initialization before the loop (useless).
I haven't run for long simulations so I'm not sure how much impact this has on the results.

--Ccarouge 15:50, 4 May 2010 (EDT)

Retrieved from "https://wiki.seas.harvard.edu/geos-chem/index.php?title=Mercury&oldid=4549"