Development of Fast-JX in GEOS-Chem

From Geos-chem
Jump to navigation Jump to search

This page describes the implementation of Fast-JX in GEOS-Chem.

A few things about Fast-JX

Jingqiu Mao implemented Fast-JX (v6.2) into GEOS-Chem v8-01-02. This version of Fast-JX includes 18 bins of wavelengths. A few good things for this new version:

  • This version can be used for stratospheric chemistry with 11 wavelength bins in UV.
  • This version of Fast-JX handles thick clouds much better, because it uses a log-spacing to add sub-layers for thick clouds or aerosol layers. We had some problems in the standard version when cloud or aerosols are too thick (the old Fast-J uses linear spacing to ad sub-layers). See Too many levels in photolysis_code.
  • This version of Fast-JX is supposed to do a better job over polar region,where solar zenith angle gets close to 90.
  • I combined most files into one module: fast_jx_mod.f. This will provide a much better platform to get updates from newer version of Fast-JX for future development.
  • Fast-JX takes into account the photolysis of O2, which will make a little bit more ozone in upper troposphere.

If you need the source code, please contact Jingqiu Mao on this.

More info on Fast-JX v6.2

Fast-JX version 6.2 corrects a long-standing problem at SZA > 89 degrees. The interpolation of the incident radiation to the mid-pt of the layer was based on the correct ray-traced solar attenuation at the edges of the layer (based on the pratmo box model grid-pt solutions). The layer edges were thus correct, but interpolating to the mid-layer was arbitrary when the layer-edge below was nearly zero (vary high SZA). The method chosen gave irregular fluctuations in the direct solar beam - and hence reported J's at mid-layer as the SZA spanned from 89 to 90+ degrees. This is now corrected with exact ray-tracing from each CTM layer mid-pt. For most SZA, differences appear at the round-off of most printed results, but near 90 degrees, errors could have been an order of magnitude, but when averaged over the day, they mostly smoothed out (from Michael Prather).

--Bob Y. 16:25, 12 May 2014 (EDT)

How does the new Fast-JX work?

Fast-JX calculates scattering at five wavelength 200nm, 300nm, 400nm, 600nm and 999nm. Fast-J calculates at four wavelength 300nm, 400nm, 600nm and 999nm (without 200nm).

  1. In chemdr.f, Fast-JX is initialized (call INPHOT(LLTROP, NPHOT)).
  2. In chemdr.f, Fast-JX is calculated by calling FAST_J( SUNCOS, OPTD, UVALBEDO).
  3. In fast_j.f, it was first determined which cloud overlap scheme (linear, random,or maximum) to be used. Then send single column properties to PHOTOJ for calculation, including:
    1. pressure profile
    2. temperature profile
    3. surface albedo
    4. Aerosol OD profile
    5. Mineral dust OD profile
    6. Cloud OD profile
  4. In calcrate.f, send J-values to smvgear by calling fjfunc.f.

Modules in fast_jx_mod.f

  • INPHOT
    • Read in labels of photolysis (call RD_JS)
    • Match Fast-JX species with Harvard species (call JV_INDEX)
    • Read in Fast-JX X-sections (spectral data) (call RD_XXX)
    • Read in T & O3 climatology (call RD_PROF)
    • Select aerosol/cloud types to be used (call SET_AER)
  • PHOTOJ
    • Set up Air, O3, BC profiles on GEOS-Chem vertical levels (call set_prof.f), this is the same as the old fast-j.
    • calculate air mass factors by calling sphere2. This is different from the old fast-j: the new AMF2 does each of the half-layers of the CTM separately, whereas the original, based on the pratmo code did the whole layers and thus calculated the ray-path to the CTM layer edges, NOT the middle.
    • add sub-layers (JXTRA) to thick cloud/aerosol layers by calling EXTRAL. This version sets up log-spaced sub-layers of increasing thickness ATAU. Given the aerosol+cloud OD/layer in visible (600 nm) calculate how to add additonal levels at top of clouds (now uses log spacing). This was linear spacing in Fast-J, so caused lots of problems for thick clouds.
    • For each wavelength bin, call OPMIE to calculate the mean actinic flux in each layer (AVGF), then save the flux into FFF (1:LPAR).
    • call JRATET to convert FFF to J-values for all species.

How does Fast-J work in GEOS-Chem v9-02 and earlier versions?

Step 1: inphot.f

Require input files (cmn_fj.h, jv_cmn.h, ratj.d, jv_spec.dat, jv_atms.dat)

  1. Error check # of layers (call ERROR_STOP)
  2. Read in labels of photolysis rates required (call RD_JS.f)
  3. translate between GEOS-Chem species nomenclature and Fast-J species nomenclature(call JV_INDEX.f)
  4. Read in JPL spectral data set (e.g. X-sections, quantum yields) (call RD_TJPL.f)
  5. Read in T & O3 climatology (cf. Nagatani/92 and McPeters/91) (call RD_PROF.f)
  6. Select Aerosol/Cloud types to be used (call SET_AER.f)

Step 2: fast_j.f

For each (NLON,NLAT) location, call subroutine PHOTOJ (in a parallel loop to compute J-values for the entire column. J-values will be stored in the common-block variable ZPJ, and will be later accessed via function FJFUNC).

  1. Read TOMS O3 columns if it's a new month(call READ_TOMSO3 from "toms_mod.f“).
  2. Cloud overlap options
  • Linear Approximation (used up to v7-04-12)->call photoj.f
  • Approximate Random Overlap (default)->call photoj.f
  • Maximum Random Overlap (MRAN, computation intensive)->call photoj.f+mmran_16.f

Step 3: photoj.f

  • Set up Air, O3, BC profiles on GEOS-Chem vertical levels (call set_prof.f)
  • Set up cloud and surface properties (Call CLDSRF.f).
  • Set up pressure levels for O3/T climatology.
  • Select appropriate monthly and latitudinal profiles
  • Apportion O3 and T on supplied climatology z* levels onto CTM levels with mass (pressure) weighting, assuming constant mixing ratio and temperature half a layer on either side of the point supplied.
  • Calculate effective altitudes using scale height at each level
  • Add aerosol column - include aerosol types here. Currently use soot water and ice; assume black carbon x-section of 10 m2/g, independent of wavelength; assume limiting temperature for ice of -40 deg C.
  • Calculate column quantities for FAST-J: monthly mean air column and monthly mean O3 column.
  • Weight the O3 column by the observed monthly mean TOMS data.
  • Interpolate O3 column to current day (TOMS data is half month time resolution). Scale monthly O3 profile to the daily O3 profile.
  • Compute actinic flux at each GEOS-CHEM vertical level (call JVALUE.f). In which, calculate air mass factors for each layer (call sphere.f). Loop over all wavelength bins to call OPMIE.f:
  • Pick nearest Mie wavelength, no interpolation.
  • For Mie code scale extinction at 1000 nm to wavelength WAVEL (QXMIE)
  • Set up total optical depth over each CTM level, DTAUX
  • Fractional extinction for Rayleigh scattering and each aerosol type
  • Define the scattering phase fn. with mix of Rayleigh(1) & Mie(MIEDX), No. of quadrature pts fixed at 4 (M__), expansion of phase fn @ 8
  • Calculate attenuated incident beam EXP(-TTAU/U0) and flux on surface
  • Take optical properties on CTM layers and convert to a photolysis level grid corresponding to layer centres and boundaries
  • Calculate cumulative total and define levels we want J-values at. Sum upwards for levels, and then downwards for Mie code readjustments.
  • SET UP FOR MIE CODE. Transpose the ascending TTAU grid to a descending ZTAU grid. This is required so that J-values can be calculated for the centre of CTM layers; the index of these layers is kept in the jndlev array.
  • Insert new levels, working downwards from the top of the atmosphere to the surface (down in 'j', up in 'k').
  • call MIESCT.f
(1) fix scattering to 4 Gauss pts = 8-stream(CALL GAUSSP.f)
(2) solve eqn of R.T. only for first-order M=1(call LEGND0.f)
(3) call BLKSLV.f to solve the block tri-diagonal system.
(4) call GEN.f and MATIN4.f.
  • Calculate J-values for all species(call JRATET.f)
  • Scale actinic flux (FFF) by Solar distance factor (SOLF)
  • With model temperature, call XSECO2.f(FLINT.f); call XSECO3.f(FLINT.f); call XSEC1D.f(FLINT.f)
  • Calculate Jvalues for O2,O3,O1D
  • Calculate remaining J-values with T-dep X-sections

Step 4: fjfunc.f

supplies J-values to SMVGEAR solver

Implementation of Fast-JX v6.2

Jingqiu Mao (Princeton) has incorporated FAST-JX v6.2 into several research versions of GEOS-Chem. This version of FAST-JX, however, has not been incorporated into the standard GEOS-Chem code. Sebastian Eastham (MIT) has added a newer version, FAST-JX v7.0, along with the UCX chemistry mechanism into GEOS-Chem v10-01.

New input files

  1. spec2008.dat: fast-J X-sections (spectral data), to replace jv_spec.dat
  2. atmos_std.dat: T & O3 climatology data, to replace jv_atms.dat
  3. chemJ2008.d: this is a tranfer map from the J's automatically calculated in fast-JX onto the names and order in the users chemistry code. This one is to replace ratj.d

Major changes

In the following sections, we describe the major changes tahtw ere made when updating FAST-J to FAST-JX:

aerosol_mod.f

old new
#include "cmn_fj.h"
#include "jv_cmn.h"
#include "fj2008.h"
RW(R) = RAA(4,IND(N)+R-1) RW(R) = RAA(5,IND(N)+R-1)
QW(R) = QAA(4,IND(N)+R-1)*FWET + QAA(4,IND(N))*(1.d0-FWET) QW(R) = QAA(5,IND(N)+R-1)*FWET + QAA(5,IND(N))*(1.d0-FWET)
WAERSL(I,J,L,N) * QAA(4,IND(N)) WAERSL(I,J,L,N) * QAA(5,IND(N))
DAERSL(I,J,L,N-1) * QAA(4,IND(N)) DAERSL(I,J,L,N-1) * QAA(5,IND(N))
REFF = 1.0D-4 * RAA(4,IND(N)) REFF = 1.0D-4 * RAA(5,IND(N))
QAA(2,IND(N)+R-1) / QAA(4,IND(N)+R-1) QAA(3,IND(N)+R-1) / QAA(5,IND(N)+R-1)

All these are because there is one more wavelength bin for mie scattering(200nm), so RAA (effective radius) and QAA(aerosol extinction coefficient) all need to be changed accordingly.

--Bob Y. 15:57, 12 May 2014 (EDT)

dust_mod.f

old new
#include "cmn_fj.h"
#include "jv_cmn.h"
#include "fj2008.h"
DUST(I,J,L,N) * QAA(4,14+N)/ ( MSDENS(N) * RAA(4,14+N) * 1.0D-6 ) DUST(I,J,L,N) * QAA(5,14+N)/ ( MSDENS(N) * RAA(5,14+N) * 1.0D-6 )
ERADIUS(JLOOP,N) = RAA(4,14+N) * 1.0D-4 ERADIUS(JLOOP,N) = RAA(5,14+N) * 1.0D-4
( ODMDUST(I,J,L,N) * QAA(2,14+N) / QAA(4,14+N) ) ( ODMDUST(I,J,L,N) * QAA(3,14+N) / QAA(5,14+N) )
DUST(I,J,L,N) * QAA(4,14+N) / ( MSDENS(N) * RAA(4,14+N) * 1.0D-6 ) DUST(I,J,L,N) * QAA(5,14+N)/ ( MSDENS(N) * RAA(5,14+N) * 1.0D-6 )
ERADIUS(JLOOP,N) = RAA(4,14+N) * 1.0D-4 ERADIUS(JLOOP,N) = RAA(5,14+N) * 1.0D-4
( ODMDUST(I,J,L,N) * QAA(2,14+N) / QAA(4,14+N) ) ( ODMDUST(I,J,L,N) * QAA(3,14+N) / QAA(5,14+N) )

--Bob Y. 15:57, 12 May 2014 (EDT)

fast_j.f, fjfunc.f, fjfunc.f

old new
USE FAST_JX_MOD, ONLY : PHOTOJ
#include "cmn_fj.h"
#include "jv_cmn.h"
#include "fj2008.h"

--Bob Y. 15:57, 12 May 2014 (EDT)

fast_jx_mod.f

This is a new file and it is the main file for Fast-JX.

--Bob Y. 15:57, 12 May 2014 (EDT)

planeflight_mod.f

old new
#include "cmn_fj.h"
#include "jv_cmn.h"
# include "fj2008.h"
S400nm = QAA(2,IND(N)+RH-1) / & QAA(4,IND(N)+RH-1) S400nm = QAA(3,IND(N)+RH-1) / & QAA(5,IND(N)+RH-1)
S400nm = QAA(2,IND(N)+RH-1) / & QAA(4,IND(N)+RH-1) S400nm = QAA(3,IND(N)+RH-1) / & QAA(5,IND(N)+RH-1)

--Bob Y. 15:57, 12 May 2014 (EDT)

Other related changes

ch3i_mod.f

add this line at the beginning;

USE FAST_JX_MOD

Because ch3i_mod.f needs fast_j to calculate something for ch3i simulations.

chemdr.f

Add this line:

USE FAST_JX_MOD,          ONLY : INPHOT

Because inphot.f is now merged into fast_jx_mod.f

--Bob Y. 16:09, 12 May 2014 (EDT)

chemistry_mod.f

Add this line:

USE FAST_JX_MOD

--Bob Y. 16:09, 12 May 2014 (EDT)

diag48_mod.f, diag49_mod.f, diag50_mod.f, diag51_mod.f

old new
USE FAST_JX_MOD, ONLY : PHOTOJ
#include "fj2008.h"

--Bob Y. 16:09, 12 May 2014 (EDT)

Files removed from v8-01-02

Filename Notes
BLKSLV.f Merged into “fastjx_mod.f”
CLDSRF.f Merged into “fastjx_mod.f”
cmn_fj.h Input files, now merged into “fj2008.h”.
EFOLD.f Removed, not even used in fast-j
FLINT.f Merged into “fastjx_mod.f” (Three-point linear interpolation function)
GAUSSP.f Merged into “fastjx_mod.f”
GEN.f Merged into “fastjx_mod.f”
inphot.f Merged into “fastjx_mod.f”
JRATET.f Merged into “fastjx_mod.f”
JVALUE.f Merged into PHOTOJ in “fastjx_mod.f”
jv_cmn.h Input files, now merged into “fj2008.h”.
jv_index.f Merged into “fastjx_mod.f”
jv_mie.h Input files, now merged into “fj2008.h”.
LEGND0.f Merged into “fastjx_mod.f”
MATIN4.f Merged into “fastjx_mod.f”
MIESCT.f Merged into “fastjx_mod.f”
NOABS.f Removed, not even used in fast-j
OPMIE.f Merged into “fastjx_mod.f”
photoj.f Merged into “fastjx_mod.f”
rd_js.f Merged into “fastjx_mod.f”
rd_prof.f Merged into “fastjx_mod.f”
RD_TJPL.f This is changed to “RD_XXX” in “fastjx_mod.f”.
set_aer.f Merged into “fastjx_mod.f”
set_prof.f Merged into “fastjx_mod.f”
SPHERE.f Merged into “fastjx_mod.f”, now changed to sphere2 function.
XSEC1D.f Merged into PHOTOJ in “fastjx_mod.f”
XSECO2.f Merged into PHOTOJ in “fastjx_mod.f”
XSECO3.f Merged into PHOTOJ in “fastjx_mod.f”

--Bob Y. 16:11, 12 May 2014 (EDT)

Some notes from Jingqiu Mao's implementation

  • The original version of Fast-JX v6.2 is based on aerosol category from UCI and U-Michigan. I have to modify accordingly to use our aerosols.
  • For our aerosols, since we don’t have data at 200 nm, I just use the data at 300 nm (which was agreed by MJP).
  • Solar distance factor, I found that this was
solf=1.d0-(0.034d0*cos(dble(iday-172)*2.d0*pi/365.d0))

now in FJX v6.2, it becomes

SOLFX  = 1.d0-(0.034d0*cos(dble(NDAY-186)*2.d0*PI/365.d0))

which one is correct?

  • The standard version use cloud OD at 600nm to determine how many sub-layers to be inserted. Here we just use the cloud optical depth from GEOS-5 met field to do that. This should be just fine. According to GEOS-5 manual, the optical depth from GEOS-5 is for the visible band (400-690 nm).
  • In old OPMIE.f, we do the scaling (For Mie code scale extinction at 1000 nm to wavelength WAVEL (QXMIE)):
      
do I=1,MX
        QXMIE(I) = QAA(KM,MIEDX(I))/QAA(4,MIEDX(I))
        SSALB(I) = SSA(KM,MIEDX(I))
enddo

Now this part is moved out of OPMIE and now into PHOTOJ, because we want to handle everything out of OPMIE, and only send DTAUX and POMGEAX to OPMIE.

  • One trick here. In aerosol_mod.f, aerosol optical depth (ODAER) is calculated at 1000 nm. When this variable is sent to Fast-J, we have to scale it up by doing this:
do I=1,MX
        QXMIE(I) = QAA(KM,MIEDX(I))/QAA(4,MIEDX(I))
        SSALB(I) = SSA(KM,MIEDX(I))
enddo

This part is not in the standard Fast-JX v6.2 (the original version). In old OPMIE.f, we then do this to scale aerosol optical depth (AER) back:

XLAER(I)=AER(I,J)*QXMIE(I)
  • In the standard Fast-JX, for scattering data, PAA(1) is always 1.0, that is why it only reads from the second term to the eighth term (only seven terms for phase function in their input file).
  • To innclude Rayleigh scattering into scattering phase function, the old fastj do this:
POMEGAX(I,J) = PIRAY(J)*PAA(I,KMIE,1)

the new fastjx do this:

POMEGAX(1,L) = POMEGAX(1,L) + 1.0d0*ODRAY
POMEGAX(3,L) = POMEGAX(3,L) + 0.5d0*ODRAY 

Note that these two treatments are actually identical, because Rayleigh phase function only has two values pi(0)=1, pi(3)=0.5.

  • In old fast-j, we commented out this line:
if( WAVE .gt. 800.d0 ) KMIE=5  ! use 999 nm prop for 800- nm

I think this is mainly for other calculations at 1000 nm in aerosol_mod.f and dust-mod.f.

  • I commented out all the code for heat flux calculations, which seem unnecessary for now.

Checklist for merging into the new version

  • acetone stuff, this was later implemented. But we can just remove that (fjx_acet_mod.f) since all of these are now in Fast-JX.
  • new optical input added by Colette and Randall? these should merged.
  • the yield of HNO4 photolysis has to be modified.

Implementation of Fast-JX v7.0

FAST-JX v7.0 has been introduced into GEOS-Chem v10-01, along with the UCX chemistry mechanism.

Overview

This update is being tested in the 1-month benchmark simulation v10-01c.

Sebastian Eastham incorporated Fast-JX v7.0a into the GEOS-Chem UCX mechanism. From Eastham et al. (2014):

GEOS-Chem uses a customized version of the Fast-JX v6.2 photolysis rate solver (Wild et al., 2000), which efficiently estimates tropospheric photolysis. The customized version uses the wavelength bands from the older Fast-J tropospheric photolysis scheme and does not consider wavelengths shorter than 289 nm, assuming they are attenuated above the tropopause. However, these high-energy photons are responsible for the release of ozone-depleting agents in the stratosphere. The standard Fast-JX model (Prather, 2012) addresses this limitation by expanding the spectrum analyzed to 18 wavelength bins covering 177–850 nm, extending the upper altitude limit to approximately 60 km. We therefore incorporate Fast-JX v7.0a into GEOS-Chem UCX. Fast-JX includes cross-section data for many species relevant to the troposphere and stratosphere. However, accurately representing sulfur requires calculation of gaseous H2SO4 photolysis, a reaction which is not present in Fast-JX but which acts as a source of sulfur dioxide in the upper stratosphere. Based on a study by Mills (2005), the mean cross-section between 412.5 and 850 nm is estimated at 2.542 × 10−25 cm2. We also add photolysis of ClOO and ClNO2, given their importance in catalytic ozone destruction, using data from JPL 10-06 (Sander et al., 2011). Fast-JX v7.0a includes a correction to calculated acetone cross sections. Accordingly, where hydroxyacetone cross-sections were previously estimated based on one branch of the acetone decomposition, a distinct set of cross sections from JPL 10-06 are used.
The base version of GEOS-Chem uses satellite observations of total ozone columns when determining ozone-related scattering and extinction. The UCX allows either this approach, as was used for the production of the results shown, or can employ calculated ozone mixing ratios instead, allowing photolysis rates to respond to changes in the stratospheric ozone layer.

New input files

  1. FJX_spec.dat: Fast-J X-sections, to replace jv_spec.dat (GEOS-Chem v9-02 and earlier versions) and spec2008.dat (Fast-JX v6.2 implementation)
  2. FJX_j2j.dat: Links GEOS-Chem chemical species to Fast-JX species, to replace ratj.d
  3. jv_spec_mie.dat: Aerosol optical properties at 5 wavelengths

Major changes

Known issues

Error in reducing wavelength bins for tropospheric chemistry only simulation

This update will be tested in the 1-month benchmark simulation v10-01c.

Sebastian Eastham wrote:

In fast_jx_mod, specifically RD_XXX, there is a transformation to reduce 18 cross sections to 12. Since bin 18 now corresponds to bin 12 and so on, the wavelengths are moved within the cross section array QQQ. However, the 12-bin capability is rarely used (if ever), so when Fast-JX was extended to allow cross sections with 1 or 3 sets of data, the 12 and 8 bin codes were not updated accordingly. This results in very large cross sections for acetone at long wavelengths, because the shorter wavelength data is being used instead.
I've notified Michael Prather - he did not know about this bug and is putting together a fix ASAP. I've written my own fix in the meantime, which results in the acetone cross sections matching much more closely, at least between the two v10-01c versions.

--Melissa Sulprizio 10:39, 12 May 2014 (EDT)

References

  1. Eastham, S.D., Weisenstein, D.K., Barrett, S.R.H., Development and evaluation of the unified tropospheric–stratospheric chemistry extension (UCX) for the global chemistry-transport model GEOS-Chem, Atmos. Env., June 2014.