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 fromenergyrelated 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)