Development of Fast-JX in 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
I 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 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).
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).
- In chemdr.f, Fast-JX is initialized (call INPHOT(LLTROP, NPHOT)).
- In chemdr.f, Fast-JX 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 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 the current Fast-J work in GEOS-Chem?
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 GEOS-Chem species nomenclature and Fast-J species nomenclature(call JV_INDEX.f)
- Read in JPL spectral data set (e.g. X-sections, 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 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).
- 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 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
New input files
- spec2008.dat: fast-J X-sections (spectral data)
- atmos_std.dat: T & O3 climatology data
- 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.
Major changes
- 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.
- 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) ) |
- fast_j.f, fjfunc.f, fjfunc.f
| USE FAST_JX_MOD, ONLY : PHOTOJ | |
| # include "cmn_fj.h" # include "jv_cmn.h" | # include "fj2008.h" |
- fast_jx_mod.f- this is a new file and it is the main file for Fast-JX.
- 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) |
- 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
- chemistry_mod.f
Add this line: USE FAST_JX_MOD
- diag48_mod.f, diag49_mod.f, diag50_mod.f, diag51_mod.f
| old | new |
|---|---|
| USE FAST_JX_MOD, ONLY : PHOTOJ | |
| # include "cmn_fj.h" # include "jv_cmn.h" | # include "fj2008.h" |
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” |
Some notes from my 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.