http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_v8-02-03

Overview 
BETA RELEASE -- October 2009 

Previous version: GEOS-Chem v8-02-02 
Next version: GEOS-Chem v8-02-04 
What's new in this version 
Contains everything in v8-02-02, plus: 

1.Kpp chemical solver package (Adrian Sandu group @ Virginia Tech) 
2.ICOADS ship emissions (Chulkyu Lee) 
3.Use of GEOS-5 ozone columns (C. Carouge) 
4.Removal of support for obsolete SGI and COMPAQ compilers 
5.Updated CH4 offline simulation (K. Wecht, C. Pickett-Heaps) 
6.Updated makefile structure (G-C support team) 
7.Correction of the time step when to calculate diagnostics (C. Carouge, P. Le Sager) 
8.Updated OTD/LIS redistribution factors for GEOS-5 
--Bob Y. 09:58, 4 February 2010 (EST) 

Previous issues now resolved in v8-02-03 
Patch for KPP for v8-02-03 
Claire Carouge wrote: 

The problem we were having with KPP in GEOS-Chem v8-02-03 is now solved. We've posted patched files on the FTP site: 
   ftp ftp.as.harvard.edu
   cd pub/geos-chem/patches/v8-02-03
   get chemdr.f
   get chemistry_mod.f
The problem is that, somehow, KPP considers 3 species declared as "dead" in globchem.dat as "active" species. As a consequence, when updating the species concentrations at the end of the chemistry, the concentrations of 3 "inactive" species were modified whereas they should stay constant. With time, this would propagate problems through out the whole chemistry. 
The fix is to update only for the number of active species as given by globchem.dat and not as given by KPP. 
The fix will be standardized in GEOS-Chem v8-02-04 and later versions. 
The reason why KPP has 3 extra active species is part of a larger discussion on the mapping of the species in KPP and was not solved yet. But the fix gives a safe way to use KPP in GEOS-Chem. 
--Bob Y. 16:54, 4 February 2010 (EST) 

Corrected Bond et al BC/OC emissions
Eric Leibensperger wrote: 

The carbon_200905 data were generated with buggy code and generated bogus files. The files in carbon_200905 SHOULD NOT be used. 
In addition, my previous statements that the Cooke overwrites over the US led to a high bias in BC and OC was WRONG. I think Cooke should still be the default (meaning in input.geos, the use Cooke line should be defaulted as T, rather than F). 
The corrected files are now contained in the carbon_200909 directories. You may obtain these from: 

ftp ftp.as.harvard.edu
cd /pub/geos-chem/data/GEOS_1x1/carbon_200909
mget *
cd /pub/geos-chem/data/GEOS_2x2.5/carbon_200909
mget *
cd /pub/geos-chem/data/GEOS_4x5/carbon_200909
mget *
GEOS-Chem v8-02-03 now points to the corrected files in carbon_200909. 

The files in the previous directory carbon_200905 have been deleted from the FTP server. The carbon_200905 directory is now a symbolic link to carbon_200909. 

Bob Yantosca has archived the carbon_200905 directories in case anyone needs to do a test with these files for backwards compability in model development. Please contact him directly if you need to use these files. It is recommended to migrate directly to the files in carbon_200909. 

Reference: 

Bond, T.C. et al.: Historical emissions of black and organic carbon aerosol fromenergy-related combustion, Global Biogeochem. Cycles, 21, GB2018, 1850-2000, doi: 10.1029/2006GB002840, 2007. 
--Bob Y. 13:32, 21 September 2009 (EDT) 

Several bug fixes in sulfate_mod.f 
Philippe Le Sager wrote: 

I realized the SO4 derived from Canadian CAC inventory was ignored in v8-02-02, because of a typo. Jenny Fisher pointed out that SO4 derived from BRAVO was not accounted for. Paul Hezel provided a fix for diagnostic ND05, which now correctly records P(SO2) from OH. Paul also found out that the P(MSA) from DMS above the PBL was not accounted for. 
All these bugs have been fixed. A patch is available for GEOS-Chem v8-02-02 and earlier versions: 
   ftp ftp.as.harvard.edu
   cd pub/geos-chem/patches/v8-02-02
   get sulfate_mod.f
This fix has been standardized in GEOS-Chem v8-02-03. 
--phs 16:57, 24 August 2009 (EDT) 

Missing NOx data in S.E.-Asia 
Data outside China from the Streets NOx inventory have been missing. The Streets_NOx_FF_2004_monthly.generic.1x1 files should be used only to provide monthly variations. 

If you are using GEOS-Chem v8-02-02 or earlier, please update your code with the patch on the ftp site: 

  ftp ftp.as.harvard.edu
  cd pub/geos-chem/patches/v8-02-02
  get streets_anthro_mod.f
This fix has been standardized in GEOS-Chem v8-02-03. 

--phs 14:14, 19 August 2009 (EDT) 

Mis-calculation of Courant numbers in tpcore_fvdas_mod.f90 
I found that the Courant numbers where not calculated for the second band of latitude and the second last one in tpcore_fvdas_mod.f90. So crap values were used instead at these locations. 

If you are using GEOS-Chem v8-02-02 or earlier, please update your code with the patch on the ftp site: 

  ftp ftp.as.harvard.edu
  cd pub/geos-chem/patches/v8-02-02
  get tpcore_fvdas_mod.f90
This fix has been standardized in GEOS-Chem v8-02-03. 

--Ccarouge 11:10, 18 August 2009 (EDT) 

Format problem in planeflight_mod.f 
Jingqiu Mao reported a problem in the plane.log output files, in which it looked like some tracers had abnormally large values. 

Philippe Le Sager wrote: 

I found a solution. The problem is that Fortran output in ASCII switches between 2 formats: 
   1.234E-12
   1.234-12
Correct, the E can disappear. And that confuses IDL. With the format specified by planeflight_mod.f (line 1237), it happens when the exponent is less than -100. 
To keep the routine general, you can force the number of digit in the exponent with a 'e' in the format descriptor as follows: 
   y = -1.233d-122
   write(6, '(es10.3)'   ), y   ! what is used now
   write(6, '(es11.3e3)' ), y   ! what we should use instead
   write(6, '(es12.3e4)' ), y
gives 
   -1.233-122
   -1.233E-122
   -1.233E-0122
So replace (es10.3) with (es11.3e3) line 1237 of planeflight_mod.f and that should fix it. 
--Bob Y. 11:35, 13 July 2009 (EDT) 

Minor fixes in wet deposition 
Bob Yantosca wrote: 

(1) In subroutine RAINOUT (in wetscav_mod.f), I think we need to apply a separate conversion factor for NH3 (i.e. as is now done in subroutine COMPUTE_F). So we'd need to change the lines: 
    REAL*8, PARAMETER    :: CONV = 8.27042925126d-1
to: 
    ! CONV_H2O2 = 0.6 * SQRT( 1.9 ), used for the ice to gas ratio for H2O2
    ! 0.6 is ( sticking  coeff H2O2  / sticking  coeff  water )
    ! 1.9 is ( molecular weight H2O2 / molecular weight water )
    REAL*8, PARAMETER :: CONV_H2O2 = 8.27042925126d-1

    ! CONV_NH3 = 0.6 * SQRT( 0.9 ), used for the ice to gas ratio for NH3
    ! 0.6 is ( sticking  coeff  NH3 / sticking  coeff  water )
    ! 0.9 is ( molecular weight NH3 / molecular weight water )
    REAL*8, PARAMETER :: CONV_NH3  = 5.69209978831d-1
and then change these acccordingly: 
    ELSE IF ( N == IDTH2O2 ) THEN

       ! Compute ice to gas ratio for H2O2 by co-condensation
       ! (Eq. 9, Jacob et al, 2000)
       IF ( C_H2O(I,J,L) > 0d0 ) THEN
          !------------------------------------------------------
          ! Prior to 7/20/09
          ! Now multiply by CONV_H2O2 (bmy, 7/20/09)
          !I2G = ( CLDICE(I,J,L) / C_H2O(I,J,L) ) * CONV
          !------------------------------------------------------
          I2G = ( CLDICE(I,J,L) / C_H2O(I,J,L) ) * CONV_H2O2
       ELSE
          I2G = 0d0
       ENDIF

    ...

    ELSE IF ( N == IDTNH3 ) THEN

       ! Compute ice to gas ratio for NH3 by co-condensation
       ! (Eq. 9, Jacob et al, 2000)
       IF ( C_H2O(I,J,L) > 0d0 ) THEN
          !--------------------------------------------------
          ! Prior to 7/20/09
          ! Now multiply by CONV_NH3 (bmy, 7/20/09)
          !I2G = ( CLDICE(I,J,L) / C_H2O(I,J,L) ) * CONV
          !--------------------------------------------------
          I2G = ( CLDICE(I,J,L) / C_H2O(I,J,L) ) * CONV_NH3
       ELSE
          I2G = 0d0
       ENDIF
Havala gave us the fix for separate CONV_H2O2 and CONV_NH3 factors for subroutine COMPUTE_F. But then we may have never thought to also apply this to subroutine RAINOUT, which subroutine, where the algorithm is similar. In any case, this is a simple fix. 
(2) In subroutine COMPUTE_F, for NH3 we have: 
          ! F is the fraction of NH3 scavenged out of the updraft
          ! (Eq. 2, Jacob et al, 2000)
          F(I,J,L) = 1d0 - EXP( -K * BXHEIGHT(I,J,L) / Vud(I,J) )
Note that BXHEIGHT(I,J,L) is used in the equation. However, almost all of the other non-aerosol tracers in COMPUTE_F use: 
          ! Distance between grid box centers [m]
          TMP = 0.5d0 * ( BXHEIGHT(I,J,L-1) + BXHEIGHT(I,J,L) )

          ! F is the fraction of NH3 scavenged out of the updraft
          ! (Eq. 2, Jacob et al, 2000)
          F(I,J,L) = 1d0 - EXP( -K * TMP / Vud(I,J) )
So in other words TMP is the distance between grid box centers rather than the distance between grid box edges. 
There is a comment: 
   !  (7 ) Now set F=0 in the first level for all tracers.  Also now
   !        compute the distance between grid box centers and use that in
   !        in Eq. 10 from Jacob et al, 2000 to compute F. (hyl, bmy, 1/24/02)
It may have been that we just never switched BXHEIGHT to TMP in the IF block for NH3 tracer, I don't recall. It is probably not terrible as-is...but we may want to change that to be consistent w/ what is used for the other tracers. 
Hongyu Liu replied: 

Regarding item (2) above, you're right that we need to be consistent throughout the code by using the distance between the grid centers. Using the BXHEIGHT makes almost no difference numerically, according to my test with radionuclide simulations conducted 10 years ago. 
--Bob Y. 10:49, 20 July 2009 (EDT) 

Minor fixes for IBM XLF compiler 
Gabriel Morin wrote: 

In ocean_mercury_mod.f, replace line 491: 
          EF     = MAX( (0.63 - 0.02 * TC), 0.0) ! keep export > 0
by the line: 
          EF     = MAX( (0.63d0 - 0.02d0 * TC), 0.0d0) ! keep export > 0
In lightning_nox_mod.f, replace line 1441: 
    CC = MAX( CCTHICK * 1d-3, 5.5 )
by the line: 
    CC = MAX( CCTHICK * 1d-3, 5.5d0 )
--Bob Y. 11:35, 13 July 2009 (EDT) 

Minor fixes in gamap_mod.f 
Dylan Millet wrote: 

There's a minor bug in gamap_mod.f. In the ND28 section, MOLC for ALD2 is given as 3 instead of 2. 
Bob Yantosca wrote: 

Junhua Liu found a bug in the gamap_mod.f. If you use GEOS-4 then the mass flux variable is "ZMMU"; however, the name "CLDMAS" gets printed out to the tracerinfo.dat file. 
These bugs have now been corrected. 

--Bob Y. 14:04, 9 October 2009 (EDT) 

Diagnostic indexing problem for timeseries output 
Amos Tai wrote: 

I just have a really quick question regarding output files... I just ran v8-02-02 with 54 tracers, saving outputs as ts_24h_avg.20050101.bpch, including ND50 diagnostics 93-99. However, as I ran GAMAP, the diaginfo/tracerinfo.dat files don't seem to give me the right variable name/unit for some variables: 
    55 :  PEDGE-$   21   PSURF 10001          hPa 175392.00(2005010400)  72 46  1
    56 : DAO-FLDS   21   PARDR 11021          hPa 175392.00(2005010400)  72 46  1
    57 : DAO-3D-$   21    UWND 12001          m/s 175392.00(2005010400)  72 46  1
    58 : DAO-3D-$   21    VWND 12002          m/s 175392.00(2005010400)  72 46  1
    59 : DAO-3D-$   21    TMPU 12003          hPa 175392.00(2005010400)  72 46  1
e.g. 56 should be SLP, but it shows PARDR; 59 should have unit of K, but it shows hPa; 63 should be OH conc, but it doesn't show the name. I'm afraid I may wrongly interpret the data. Do you have a quick answer on what's wrong? 
Bob Yantosca wrote: 

Philippe & I looked into this. We discovered a discrepancy in the tracer #'s in the diag50_mod.f under the ND67 diagnostic category. Evidently we updated the tracer #'s in gamap_mod.f but not in the diag50_mod.f. The correct diagnostic #'s are: 
   HFLUX    GMAO HFLUX field               0.000E+00  1    11001 1.000E+00 W/m2
   RADSWG   GMAO RADSWG field              0.000E+00  1    11002 1.000E+00 W/m2
   PREACC   GMAO PREACC field              0.000E+00  1    11003 1.000E+00 W/m2
   PRECON   GMAO PRECON field              0.000E+00  1    11004 1.000E+00 W/m2
   TS       GMAO TS field                  0.000E+00  1    11005 1.000E+00 K
   RADSWT   GMAO RADSWT field              0.000E+00  1    11006 1.000E+00 W/m2
   USTAR    GMAO USTAR field               0.000E+00  1    11007 1.000E+00 m/s
   Z0       GMAO Z0 field                  0.000E+00  1    11008 1.000E+00 m
   PBL      GMAO PBL field                 0.000E+00  1    11009 1.000E+00 hPa
   CLDFRC   GMAO CLDFRC field              0.000E+00  1    11010 1.000E+00 unitless
   U10M     GMAO U10M field                0.000E+00  1    11011 1.000E+00 m/s  
   V10M     GMAO V10M field                0.000E+00  1    11012 1.000E+00 m/s
   PS-PBL   GMAO PS-PBL field              0.000E+00  1    11013 1.000E+00 hPa
   ALBD     GMAO ALBD field                0.000E+00  1    11014 1.000E+00 unitless
   PHIS     GMAO PHIS field                0.000E+00  1    11015 1.000E+00 m
   CLDTOP   GMAO CLDTOP field              0.000E+00  1    11016 1.000E+00 level
   TROPP    GMAO TROPP field               0.000E+00  1    11017 1.000E+00 hPa
   SLP      GMAO SLP field                 0.000E+00  1    11018 1.000E+00 hPa
   TSKIN    GMAO TSKIN field               0.000E+00  1    11019 1.000E+00 K
   PARDF    GMAO PARDF field               0.000E+00  1    11020 1.000E+00 W/m2
   PARDR    GMAO PARDR field               0.000E+00  1    11021 1.000E+00 W/m2
   GWET     GMAO GWET field                0.000E+00  1    11022 1.000E+00 unitless
   EFLUX    GMAO EFLUX field               0.000E+00  1    11023 1.000E+00 W/m2
But in diag50_mod.f we have: 
       ELSE IF ( N == 95 ) THEN

          !---------------------
          ! Sea level prs [hPa]
          !---------------------                      CATEGORY = 'DAO-FLDS'
          UNIT     = 'hPa'
          GMNL     = 1
          GMTRC    = 21
so the tracer # for SLP in diag50_mod.f is being written out as #21 but it should be #18. GAMAP will use the tracer # and the units from the bpch file if they are present (if not, it refers to diaginfo.dat & tracerinfo.dat). That is the problem. 
Long story short: historical baggage. 
What you can do for the time being is to manually edit the line in 
   ~tai/runs/run.v8-02-02.20091009/tracerinfo.dat
that says: 
   PARDR    GMAO PARDR field               0.000E+00  1    11021 1.000E+00 W/m2
to 
   #PARDR    GMAO PARDR field               0.000E+00  1    11021 1.000E+00 W/m2
   SLP      GMAO SLP field                 0.000E+00  1    11021 1.000E+00 hPa
(The # is the comment character.) This will make it pick up the #21 in the bpch file w/o you having to re-run the code. 
Also, in diag50_mod.f, the unit string was not defined, which caused the units of the temperature field to be displayed incorrectly. The following fix was made: 
        ELSE IF ( N == 99 ) THEN

           !---------------------
           ! Temperature
           !---------------------
           CATEGORY  = 'DAO-3D-$'
           UNIT      = 'K'          ! Add unit string (bmy, 10/13/09)
           GMNL      = ND50_NL
           GMTRC     = 3
This has now been fixed in GEOS-Chem v8-02-03. 
--Bob Y. 10:07, 13 October 2009 (EDT) 

Avoiding the "Too many levels in photolysis code" error 
You will recall the following error message from the FAST-J photolysis code. In order to avoid this error, we have increased the NL parameter from 750 to 1000 in header file jv_mie.h: 

!-----------------------------------------------------------------------
!     NL=750 was too small again, so we upped it to 1000. 
!     Uncomment this line to restore the previous definition (phs, 10/9/09)
!     PARAMETER (NL=750, N__=2*NL, M__=4)
!-----------------------------------------------------------------------
     PARAMETER (NL=1000, N__=2*NL, M__=4)
In some instances, the Gaussian quadrature algorithm of FAST-J may require more than 750 points. This is not necessarily an error, but can be a reflection of local conditions. Several users have reported encountering this error recently, so we decided to increase the NL parameter to try to avoid the error. 

The error only is an issue because FAST-J was written with old-style F77 COMMON blocks, which require a maximum size. Long story short: historical baggage. 

--Bob Y. 10:15, 13 October 2009 (EDT) 

Outstanding issues not yet resolved in v8-02-03 
Memory issue with KPP 
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. 

In the mean time, if you experience a memory problem with GEOS-Chem version 8-02-03 with KPP (resulting in unexpected crash of the run), please revert to SMVGEAR. For this, you need to change one line in the chemistry menu in input.geos: 

USE solver coded by KPP : T
replaced by: 

USE solver coded by KPP : F
--Ccarouge 12:23, 5 November 2009 (EST) 

KPP is not compatible with IFORT 9.1 
Aaron van Donkelaar wrote: 

I've just downloaded GEOS-Chem v8-02-03 and can't get it to run with the KPP solver (SMVGEAR appears to run fine). It crashes after several errors similar to: 
-----------------------------------------------------

 Forced exit from Rosenbrock due to the following error:
 --> Step size too small: T + 10*H = T or H < Roundoff
 T=  0.000000000000000E+000 and H=  1.000000000000000E-005
 failed twice !!!
 
No. of function calls:     0
No. of jacobian calls:     0
No. of steps:              0
No. of accepted steps:     0
No. of rejected steps      0
       (except at very beginning)
No. of LU decompositions:                  0
No. of forward/backward substitutions:     0
No. of singular matrix decompositions:     0
 
Texit, the time corresponding to the
       computed Y upon return:                 0.0000
Hexit, last accepted step before exit:         0.0000
Hnew, last predicted step (not yet taken):     0.0000
 
 JLOOP, I, J, L            5           5           1           1
 Forced exit from Rosenbrock due to the following error:
 --> Step size too small: T + 10*H = T or H < Roundoff
 T=  0.000000000000000E+000 and H=  1.000000000000000E-005
 failed twice !!!
-------------------------------------------------------
 
I recompiled the solver, but the result is the same. Any advice? 
Bob Yantosca wrote: 

I confirmed the same error you got with IFORT 9.1. Probably the best thing to do is to move to IFORT 10.1. The code runs without error on IFORT 10.1. We were able to compile and run with this compiler version: 
   Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 10.1
   Build 20080212 Package ID: l_fc_p_10.1.013
   Copyright (C) 1985-2008 Intel Corporation.  All rights reserved.
It's probably not worth trying to fix this, as IFORT 9.1 is now an ancient compiler version. Also, IFORT 10 gives you better performance and optimization anyway. 
Aaron van Donkelaar replied: 

KPP works fine with IFORT 10.0.023. Seems to have been a compiler issue. 
IFORT 11 
UPDATE: As of March 2010, GEOS-Chem is now compatible with Intel Fortran Compiler (aka IFORT) version 11.1.058 or higher versions. Bob Yantosca has done some timing tests with IFORT version 11.1.069. 

--Bob Y. 16:36, 25 March 2010 (EDT) 

ND65 diagnostic and KPP solver 
Be advised, the ND65 prod/loss diagnostic will probably not work correctly when using the KPP solver. This is because the ND65 diagnostic arrays were manually written into the smvgear.f routine. Becuase the KPP solver code is automatically generated by a script (supplied to us by the Virginia Tech group), it was not possible to insert the ND65 arrays into the KPP routines. 

Therefore, if you wish to obtain the ND65 prod/loss output, you should continue to use the SMVGEAR solver in GEOS-Chem. 

XLF compiler not yet tested 
The GEOS-Chem support team has not tested the new GEOS-Chem Makefile Structure using the IBM XLF compiler. Users who run GEOS-Chem on IBM machines are asked to inform the GEOS-Chem support team of any bugs or issues. 

--Bob Y. 16:39, 13 October 2009 (EDT) 

