GEOS-Chem Newsletter December 2009
Dear GEOS-Chem Users,
Here is your December 2009 GEOS-Chem Newsletter. We would like to take this opportunity to wish everyone Happy Holidays!
Thank you for your continued support of GEOS-Chem,
Bob Yantosca and Claire Carouge
- 1 GEOS-Chem model management and logistics
- 2 GEOS-Chem Support Team news
- 3 GEOS-Chem update
- 4 Column code update
GEOS-Chem model management and logistics
GEOS-Chem Steering Committee
The last GEOS-Chem Steering Committee telecon took place on December 8, 2009. Minutes of this meeting (and past meetings) are available at HERE.
The next Steering Committee Telecon shall take place on Tuesday, March 9, 2010 at 10:30 EST. Mark your calendars!
GEOS-Chem Working Groups
Please see the GEOS-Chem Working Groups web page for links to each of the Working Groups.
We encourage each GEOS-Chem user to join the Working Group that is most relevant to his/her area of research.
Please take a moment to review the GEOS-Chem website. Please inform the GEOS-Chem Support Team of any new developments such as newly submitted or published papers, personnel changes in your user group, updated research blurbs, etc.
ConfirmAccount has been installed on the GEOS-Chem wiki. This allows the GEOS-Chem support team to manually confirm all requests for new wiki accounts. This is a measure intended to defeat spamming and unauthorized modification of wiki pages.
GEOS-Chem Support Team news
With Philippe's departure last month, the GEOS-Chem Support Team now consists of Bob Yantosca and Claire Carouge. We are working hard to bring many new updates into GEOS-Chem v8-02-04 and GEOS-Chem v8-02-05!
Please note that due to circumstances, there may be some times when we may not be able to respond to your queries right away. We ask for your patience...we will attempt to answer your questions as soon as possible!
We would like to ask all GEOS-Chem users to be resourceful, especially when it comes to coding and debugging. Please see our GEOS-Chem Style Guide for some helpful information about coding practices.
Coding and debugging Tips
We recommend that all users practice defensive programming techniques; that is, to write code in such a way to avoid some of the most common pitfalls, such as:
Division by zero
One error that we see quite frequently is the division-by-zero error. This often occurs when code that was originally written for one context is modified or extended to be used in a different context.
For example, in GEOS-Chem we recently found a bug of this type in source code routine arsl1k.f, which is called from calcrate.f. Function arsl1k.f has a division such as:
! Compute ARSL1K according to the formula listed above ARSL1K = ... 2.749064E-4*SQM/(STKCF*STK)
where the variable STKCF is in a denominator. A later modification to GEOS-Chem in the calling routine calcrate.f caused the value of STKCF = 0 to be passed to function arsl1k.f under certain circumstances. Thus, the same code that used to work correctly under all conditions now results in a division-by-zero error under some conditions.
Therefore, whenever you are writing code with divisions, you should ask yourself "Is there any possibility that the denominator could be set to zero, either in this routine, or in the calling routine?" If the answer is yes, then you can take the following steps:
(1) Manually set up an IF statement to bracket the division. In the above example, the following code is sufficient to avoid the error:
IF ( ... STCKF < 1d-30 ... ) THEN
ARSL1K = some alternate default value ELSE ARSL1K = ... 2.749064E-4*SQM/(STKCF*STK) ENDIF
In the GEOS-Chem code, STCKF is typically of order 0.1. If STCKF ever gets less than the "epsilon" value of 1d-30, then it is essentially zero and should not be used in the division.
This approach works well if you are reasonably certain of the order of magnitude of the denominator. Otherwise, you can take one of the following approaches:
(2) You can also use the function IS_SAFE_DIV in error_mod.f to test if a division is safe to be performed. IS_SAFE_DIV tests the value of the numerator and denominator to make sure that the quotient can be a representable floating-point number. Here is an example:
USE ERROR_MOD, ONLY : IS_SAFE_DIV ... ! Only do division if it's safe IF ( IS_SAFE_DIV( NUM, DEN ) ) THEN QUOT = NUM / DEN ENDIF
(3) You can also use the function SAFE_DIV, also in error_mod.f to perform the division. If the division cannot be performed, an user-defined alternate value is returned. For example:
USE ERROR_MOD, ONLY : SAFE_DIV ... ! Do the division. If division isn't possible ! then return the alternate value of -1 QUOT = SAFE_DIV( NUM, DEN, -1d0 ) ! Print a message if QUOT is set to the alternate value IF ( QUOT < 0d0 ) THEN PRINT*, '### The division was not performed' ENDIF
NOTE: The alternate return value largely depends on the context of the code. For example, if the quotient is always supposed to be a positive number, then an alternate value of -1 can be used to unambiguously denote that the division could not be performed.
Whether you use IS_SAFE_DIV or SAFE_DIV is largely a matter of personal preference.
For more information, please see our Floating point math issues wiki page.
Parallel DO-loop errors
Another error we see frequently is when variables are added to a parallel DO loop but are not held PRIVATE. Let's say we start with this example:
!$OMP PARALLEL DO !$OMP+SHARED( DEFAULT ) !$OMP+PRIVATE( I, J, L, N ) DO N = 1, N_TRACERS DO L = 1, LLPAR DO J = 1, JJPAR DO I = 1, IIPAR ! Update tracers STT(I,J,L,N) = ... ENDDO ENDDO ENDDO ENDDO !$OMP END PARALLEL DO
where the tracer array STT(I,J,L,N) is assigned to something. In this case STT does not need to be held PRIVATE because it has the same scope as the loop we are parallelizing over. In other words, we are parallelizing over (I,J,L,N) and STT is dimensioned for (I,J,L,N).
However, let's say we wanted to add a scalar value to this loop, such as:
!$OMP PARALLEL DO !$OMP+SHARED( DEFAULT ) !$OMP+PRIVATE( I, J, L, N ) DO N = 1, N_TRACERS DO L = 1, LLPAR DO J = 1, JJPAR DO I = 1, IIPAR ! Update tracers STT(I,J,L,N) = ...
! Convert units from [kg] to mixing ratio [v/v] TMP = STT(I,J,L,N) * TCVV(N) / AD(I,J,L)
ENDDO ENDDO ENDDO ENDDO !$OMP END PARALLEL DO
Here, the variable TMP needs to be held PRIVATE, so that each processor will get its own local copy of it. This will prevent TMP from being overwritten with junk values by the other procesors.
At the very least, our computation will be incorrect. However, we could also encounter a segmentation fault if TMP is set to zero or a very low value (i.e. 1e-315) and then ends up in an exponential, logarithm, or denominator of a division later on in the loop.
The fix is to simply add TMP to the PRIVATE statement:
!$OMP+PRIVATE( I, J, L, N, TMP )
Long story short: before adding variables into a parallel loop, you need to decide if they have to be held PRIVATE. For more information on this topic please see our Parallelizing GEOS-Chem wiki page.
The GEOS-Chem Support Team will be away from the office during the Winter Holiday Break (i.e. week between Christmas & New Year's). We shall return to the office on Monday, January 4, 2010. Please plan accordingly.
Status of GEOS-Chem v8-02-04
GEOS-Chem v8-02-04 "beta" is currently undergoing both undergoing 1-month and 1-year benchmark test simulations.
We shall run 3 1-year benchmarks:
- Run 0: "Baseline" run (with EPA/NEI05 emissions)
- Run 1: Same as Run #1, but with Linoz strat Ox chemistry turned on
- Run 2: Same as Run #2, but with updated MEGAN emissions
As of 12/17/09, Due to some technical issues we have had to restart Run0. At this point we anticipate being finished with the benchmarking process shortly after the holidays.
Also, for the first time, the output from the 1-yr benchmark simulations will contain plots of surface concentrations and profiles of PM2.5 from the IMPROVE network of monitoring sites. We would like to thank Colette Heald and the Aerosols Working Group for supplying these data.
KPP memory issue is now fixed
As you may recall from GEOS-Chem v8-02-03, running a GEOS-Chem 2 x 2.5 simulation with KPP could potentially take more than 8 GB of memory.
- It appears that KPP solver is asking for more memory than SMVGEAR. A run with KPP would ask almost twice the memory of a same run with SMVGEAR. We are working on this memory issue.
The interface between GEOS-Chem and KPP was modified in order to correct this issue. The memory load of a GEOS-Chem simulation with KPP is now only slightly higher than with SMVGEAR. The results are not affected at all by the change.
In GEOS-Chem v8-02-04 you now have the option of using the new EPA/NEI 2005 North American emissions inventory instead of the default EPA/NEI 1999 inventory.
Aaron van Donkelaar and Philippe Le Sager worked on implementing this into GEOS-Chem. For more detailed information about this inventory, please visit our EPA/NEI05 North American emissions wiki page.
Parallel make capability
GEOS-Chem v8-02-04 can now take advantage of GNU Make's "parallel make" option. This means that you can compile GEOS-Chem on more than one processor, for the purpose of reducing overall computation time.
The parallel compilation option is specified with the -j switch:
where N is the number of processors. For example, if you have 4 available processors, you can type:
then the GNU Make utility will compile 4 files at a time. When one file finishes compiling, the next available processor will grab the next file to be compiled.
The amount of speedup in compilation will depend on your system. Here are the results of some timing tests:
|Test by||Compiler||NCPU||Wall Time (mm:ss)||% speedup|
NOTE: The number of processors you use to compile GEOS-Chem is independent of the number of processors that you will use to run GEOS-Chem (which is usually set by your queue or by the environment variable OMP_NUM_THREADS).
In the pipeline
The following updates are slated to be incorporated into the GEOS-Chem v8-02-05 release.
|TOMAS aerosol Microphysics||Win T.
|Final shakedown testing & benchmarking is now being done @ CMU.|
|Updated Isoprene Chemistry||Fabien Poulot
P. Le Sager
|In beta testing. The GEOS-Chem KPP source code is currently being recompiled for the new isoprene mechanism.|
|ISORROPIA II||Havala Pye
|Code has been delivered to G-C support team|
|Updated aerosol optics||Randall Martin, Colette Heald||Updates have been received by the G-C support team.|
|Modification to SOA formation||Aerosols Working Group||Primary organic aerosol concentrations [ug/m3]
SOA formation should not depend on the mass on inorganics. Specifically, in carbon_mod.f, MPOC in the code should only consist of organic aerosol.
We carry carbon mass only in the STT array and here multiply by 2.1 to account for non-carbon mass in the SOA.
Partitioning theory (Pankow, 1994) describes organic phase partitioning assuming absorption into pre-existing organic mass. There is currently no theoretical or laboratory support for absorption of organics into inorganics.
Note that previous versions of the standard code (v7-04-07 through v8-02-04) did include absorption into inorganics. (Colette Heald, 12/3/09)</tt>
For a complete list of outstanding code updates, please see our GEOS-Chem model development priorities page.
Column code update
The GEOS-Chem column code continues to progress quickly. Bob Yantosca and Arlindo da Silva from GSFC are working on interfacing the GEOS-Chem column code to the GEOS-5 GCM.
Also, Harvard and GSFC will collaborate on the construction of a common emissions component that can be shared by both the columnized versions the GEOS-Chem and GMI models. This will also allow easy interchange of GEOS-Chem emissions into GMI and vice-versa.
We have listed our guidelines for writing columnized code on the wiki. These are good rules of thumb for any programming project, and we encourage all GEOS-Chem users to follow these guidelines whenever possible.
You can track our progress on the GEOS-Chem column code wiki page. Check back often for the latest developments!
--Bob Y. 14:51, 17 December 2009 (EST)