FASTJX v6.2 photolysis mechanism
This page describes the implementation of the FastJX photolysis code into GEOSChem. There have beeen two lines of FastJX development:
 FASTJX v6.2 has been incorporated into several research version of GEOSChem by Jingqiu Mao (Princeton)
 FASTJX v7.0 has been incorporated into GEOSChem v1001 (concurrent with the UCX chemistry mechanism) by Sebastian Eastham (MIT).
In the following sections we shall describe the development of FASTJX v6.2.
Contents
 1 A few things about FastJX v6.2
 2 How does FastJX v6.2 work?
 3 How does FastJ work in GEOSChem v902 and earlier versions?
 4 Implementation of FastJX v6.2
A few things about FastJX v6.2
Jingqiu Mao implemented FastJX (v6.2) into GEOSChem v80102. This version of FastJX 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 FastJX handles thick clouds much better, because it uses a logspacing to add sublayers for thick clouds or aerosol layers. We had some problems in the standard version when cloud or aerosols are too thick (the old FastJ uses linear spacing to ad sublayers). See Too many levels in photolysis_code.
 This version of FastJX 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 FastJX for future development.
 FastJX 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 FastJX v6.2
FastJX version 6.2 corrects a longstanding problem at SZA > 89 degrees. The interpolation of the incident radiation to the midpt of the layer was based on the correct raytraced solar attenuation at the edges of the layer (based on the pratmo box model gridpt solutions). The layer edges were thus correct, but interpolating to the midlayer was arbitrary when the layeredge below was nearly zero (vary high SZA). The method chosen gave irregular fluctuations in the direct solar beam  and hence reported J's at midlayer as the SZA spanned from 89 to 90+ degrees. This is now corrected with exact raytracing from each CTM layer midpt. For most SZA, differences appear at the roundoff 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 FastJX v6.2 work?
FastJX calculates scattering at five wavelength 200nm, 300nm, 400nm, 600nm and 999nm. FastJ calculates at four wavelength 300nm, 400nm, 600nm and 999nm (without 200nm).
 In chemdr.f, FastJX is initialized (call INPHOT(LLTROP, NPHOT)).
 In chemdr.f, FastJX is calculated by calling FAST_J( SUNCOS, OPTD, UVALBEDO).
 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:
 pressure profile
 temperature profile
 surface albedo
 Aerosol OD profile
 Mineral dust OD profile
 Cloud OD profile
 In calcrate.f, send Jvalues to smvgear by calling fjfunc.f.
Modules in fast_jx_mod.f
 INPHOT
 Read in labels of photolysis (call RD_JS)
 Match FastJX species with Harvard species (call JV_INDEX)
 Read in FastJX Xsections (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 GEOSChem vertical levels (call set_prof.f), this is the same as the old fastj.
 calculate air mass factors by calling sphere2. This is different from the old fastj: the new AMF2 does each of the halflayers of the CTM separately, whereas the original, based on the pratmo code did the whole layers and thus calculated the raypath to the CTM layer edges, NOT the middle.
 add sublayers (JXTRA) to thick cloud/aerosol layers by calling EXTRAL. This version sets up logspaced sublayers 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 FastJ, 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 Jvalues for all species.
How does FastJ work in GEOSChem v902 and earlier versions?
Step 1: inphot.f
Require input files (cmn_fj.h, jv_cmn.h, ratj.d, jv_spec.dat, jv_atms.dat)
 Error check # of layers (call ERROR_STOP)
 Read in labels of photolysis rates required (call RD_JS.f)
 translate between GEOSChem species nomenclature and FastJ species nomenclature(call JV_INDEX.f)
 Read in JPL spectral data set (e.g. Xsections, quantum yields) (call RD_TJPL.f)
 Read in T & O3 climatology (cf. Nagatani/92 and McPeters/91) (call RD_PROF.f)
 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 Jvalues for the entire column. Jvalues will be stored in the commonblock variable ZPJ, and will be later accessed via function FJFUNC).
 Read TOMS O3 columns if it's a new month(call READ_TOMSO3 from "toms_mod.f“).
 Cloud overlap options
 Linear Approximation (used up to v70412)>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 GEOSChem 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 xsection of 10 m2/g, independent of wavelength; assume limiting temperature for ice of 40 deg C.
 Calculate column quantities for FASTJ: 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 GEOSCHEM 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 Jvalues 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 Jvalues 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 = 8stream(CALL GAUSSP.f)
 (2) solve eqn of R.T. only for firstorder M=1(call LEGND0.f)
 (3) call BLKSLV.f to solve the block tridiagonal system.
 (4) call GEN.f and MATIN4.f.
 Calculate Jvalues 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 Jvalues with Tdep Xsections
Step 4: fjfunc.f
supplies Jvalues to SMVGEAR solver
Implementation of FastJX v6.2
Jingqiu Mao (Princeton) has incorporated FastJX v6.2 into several research versions of GEOSChem. This version of FastJX, however, has not been incorporated into the standard GEOSChem code. Sebastian Eastham (MIT) has added a newer version, FastJX v7.0 (along with the UCX chemistry mechanism) into GEOSChem v1001.
Input files for FASTJ v6.2
FASTJ v6.2 uses the following input files:
File  Description 

spec2008.dat 

atmos_std.dat 

chemJ2008.d 

Bob Y. 13:48, 20 May 2014 (EDT)
Major changes
In the following sections, we describe the major changes tahtw ere made when updating FASTJ to FASTJX:
aerosol_mod.f
old  new 

#include "cmn_fj.h" #include "jv_cmn.h" 
#include "fj2008.h" 
RW(R) = RAA(4,IND(N)+R1)  RW(R) = RAA(5,IND(N)+R1) 
QW(R) = QAA(4,IND(N)+R1)*FWET + QAA(4,IND(N))*(1.d0FWET)  QW(R) = QAA(5,IND(N)+R1)*FWET + QAA(5,IND(N))*(1.d0FWET) 
WAERSL(I,J,L,N) * QAA(4,IND(N))  WAERSL(I,J,L,N) * QAA(5,IND(N)) 
DAERSL(I,J,L,N1) * QAA(4,IND(N))  DAERSL(I,J,L,N1) * QAA(5,IND(N)) 
REFF = 1.0D4 * RAA(4,IND(N))  REFF = 1.0D4 * RAA(5,IND(N)) 
QAA(2,IND(N)+R1) / QAA(4,IND(N)+R1)  QAA(3,IND(N)+R1) / QAA(5,IND(N)+R1) 
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.0D6 )  DUST(I,J,L,N) * QAA(5,14+N)/ ( MSDENS(N) * RAA(5,14+N) * 1.0D6 ) 
ERADIUS(JLOOP,N) = RAA(4,14+N) * 1.0D4  ERADIUS(JLOOP,N) = RAA(5,14+N) * 1.0D4 
( 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.0D6 )  DUST(I,J,L,N) * QAA(5,14+N)/ ( MSDENS(N) * RAA(5,14+N) * 1.0D6 ) 
ERADIUS(JLOOP,N) = RAA(4,14+N) * 1.0D4  ERADIUS(JLOOP,N) = RAA(5,14+N) * 1.0D4 
( 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 FastJX.
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)+RH1) / & QAA(4,IND(N)+RH1)  S400nm = QAA(3,IND(N)+RH1) / & QAA(5,IND(N)+RH1) 
S400nm = QAA(2,IND(N)+RH1) / & QAA(4,IND(N)+RH1)  S400nm = QAA(3,IND(N)+RH1) / & QAA(5,IND(N)+RH1) 
Bob Y. 15:57, 12 May 2014 (EDT)
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 v80102
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 fastj 
FLINT.f  Merged into “fastjx_mod.f” (Threepoint 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 fastj 
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 FastJX v6.2 is based on aerosol category from UCI and UMichigan. 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(iday172)*2.d0*pi/365.d0))
now in FJX v6.2, it becomes
SOLFX = 1.d0(0.034d0*cos(dble(NDAY186)*2.d0*PI/365.d0))
which one is correct?
 The standard version use cloud OD at 600nm to determine how many sublayers to be inserted. Here we just use the cloud optical depth from GEOS5 met field to do that. This should be just fine. According to GEOS5 manual, the optical depth from GEOS5 is for the visible band (400690 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 FastJ, 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 FastJX 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 FastJX, 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 fastj, 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 dustmod.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 FastJX.
 new optical input added by Colette and Randall? these should merged.
 the yield of HNO4 photolysis has to be modified.