Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Planeflight AOD diagnostics appear to be malfunctioning (GC 12.7) [BUG/ISSUE] #294

Closed
JFBrewer opened this issue May 5, 2020 · 15 comments
Closed
Assignees
Labels
category: Bug Something isn't working
Milestone

Comments

@JFBrewer
Copy link

JFBrewer commented May 5, 2020

Hi folks,

I've been playing around with the planeflight diagnostic outputs recently, and I think I've found a bug. On this wiki page, it suggests that the planeflight diagnostic can output 6 different AOD variables (dust, sulfate, black carbon, organic carbon, sea salt in 2 modes) for both the column and below the aircraft. However, when I output these variables, all 12 options give the exact same value at every point. It would appear to me that they may not be being properly updated by the model. I'm not sure whether the given value is a total AOD value or if it's just one of the 12 diagnostics being applied in all cases.

I've attached one of my Planeflight.dat files and the corresponding output file from May 1, 2016. This was created using GC v12.7.1 and a custom quarter degree nested run set up over Korea.

Planeflight_Bug_Report.zip

Jared

@JFBrewer JFBrewer added the category: Bug Something isn't working label May 5, 2020
@yantosca
Copy link
Contributor

yantosca commented May 6, 2020

Thanks for writing. We will look into this. I wonder if this is somehow related to #293.

@yantosca yantosca self-assigned this May 6, 2020
@yantosca
Copy link
Contributor

yantosca commented May 8, 2020

I took a look at the planeflight code (see this permalink).

!--------------------------
! Column aerosol optical
! depths [unitless]
!--------------------------
CASE ( 2001:2006 )
! Remove MISSING flag
VARI(V) = 0e+0_fp
! Aerosol number
N = PVAR(V) - 2000
! Loop over RH bins
DO ISPC= 1, NAER
IF ( .not. LINTERP ) THEN
DO LL = 1, State_Grid%NZ
! Accumulate
VARI(V) = VARI(V) + ODAER(I,J,LL,IWVSELECT(1,1),ISPC)
ENDDO
ELSE
DO LL = 1, State_Grid%NZ
! Interpolated using angstrom exponent between
! Closest available wavelengths
! (coefs pre-calculated in CALC_AOD)
!catch any zero values before interpolation
IF ((ODAER(I,J,LL,IWVSELECT(2,1),ISPC).GT.0).AND. &
(ODAER(I,J,LL,IWVSELECT(1,1),ISPC).GT.0)) THEN
VARI(V) = VARI(V) + &
(ODAER(I,J,LL,IWVSELECT(2,1),ISPC)*ACOEF_WV(1)** &
(BCOEF_WV(1)*LOG(ODAER(I,J,LL,IWVSELECT(1,1),ISPC)/&
ODAER(I,J,LL,IWVSELECT(2,1),ISPC))))
ENDIF
ENDDO
ENDIF
ENDDO
!now add in the dust
IF ( .not. LINTERP ) THEN
DO LL = 1, State_Grid%NZ
DO ISPC = 1, NDUST
! Accumulate
VARI(V) = VARI(V) + ODMDUST(I,J,LL,IWVSELECT(1,1),ISPC)
ENDDO
ENDDO
ELSE
DO LL = 1, State_Grid%NZ
! Interpolated using angstrom exponent between
! Closest available wavelengths
! (coefs pre-calculated in CALC_AOD)
!catch any zero values before interpolation
DO ISPC = 1, NDUST
IF ((ODAER(I,J,LL,IWVSELECT(2,1),ISPC).GT.0).AND. &
(ODAER(I,J,LL,IWVSELECT(1,1),ISPC).GT.0)) THEN
VARI(V) = VARI(V) + &
(ODMDUST(I,J,LL,IWVSELECT(2,1),ISPC)*ACOEF_WV(1)** &
(BCOEF_WV(1)*LOG(ODMDUST(I,J,LL,IWVSELECT(1,1),ISPC)/&
ODMDUST(I,J,LL,IWVSELECT(2,1),ISPC))))
ENDIF
ENDDO
ENDDO
ENDIF
!--------------------------
! Aerosol optical depths
! below plane [unitless]
!--------------------------
CASE ( 3001:3006 )
! Remove MISSING flag
VARI(V) = 0e+0_fp
! Aerosol number
N = PVAR(V) - 3000
! Loop over RH bins
DO ISPC= 1, NAER
IF ( .not. LINTERP ) THEN
DO LL = 1, State_Grid%NZ
! Skip non-tropospheric boxes
IF ( .not. State_Met%InTroposphere(I,J,L) ) CYCLE
! Accumulate
VARI(V) = VARI(V) + ODAER(I,J,LL,IWVSELECT(1,1),ISPC)
ENDDO
ELSE
DO LL = 1, State_Grid%NZ
! Skip non-tropospheric boxes
IF ( .not. State_Met%InTroposphere(I,J,L) ) CYCLE
! Interpolated using angstrom exponent between
! Closest available wavelengths
! (coefs pre-calculated in CALC_AOD)
!catch any zero values before interpolation
IF ((ODAER(I,J,LL,IWVSELECT(2,1),ISPC).GT.0).AND. &
(ODAER(I,J,LL,IWVSELECT(1,1),ISPC).GT.0)) THEN
VARI(V) = VARI(V) + &
(ODAER(I,J,LL,IWVSELECT(2,1),ISPC)*ACOEF_WV(1)** &
(BCOEF_WV(1)*LOG(ODAER(I,J,LL,IWVSELECT(1,1),ISPC)/&
ODAER(I,J,LL,IWVSELECT(2,1),ISPC))))
ENDIF
ENDDO
ENDIF
ENDDO
!now add in the dust
IF ( .not. LINTERP ) THEN
DO LL = 1, State_Grid%NZ
! Skip non-tropospheric boxes
IF ( .not. State_Met%InTroposphere(I,J,L) ) CYCLE
DO ISPC = 1, NDUST
! Accumulate
VARI(V) = VARI(V) + ODMDUST(I,J,LL,IWVSELECT(1,1),ISPC)
ENDDO
ENDDO
ELSE
DO LL = 1, State_grid%NZ
! Skip non-tropospheric boxes
IF ( .not. State_Met%InTroposphere(I,J,L) ) CYCLE
! Interpolated using angstrom exponent between
! Closest available wavelengths
! (coefs pre-calculated in CALC_AOD)
!catch any zero values before interpolation
DO ISPC = 1, NDUST
IF ((ODAER(I,J,LL,IWVSELECT(2,1),ISPC).GT.0).AND. &
(ODAER(I,J,LL,IWVSELECT(1,1),ISPC).GT.0)) THEN
VARI(V) = VARI(V) + &
(ODMDUST(I,J,LL,IWVSELECT(2,1),ISPC)*ACOEF_WV(1)** &
(BCOEF_WV(1)*LOG(ODMDUST(I,J,LL,IWVSELECT(1,1),ISPC)/&
ODMDUST(I,J,LL,IWVSELECT(2,1),ISPC))))
ENDIF
ENDDO
ENDDO
ENDIF

I think there are a few things wrong here:

First, the aerosols & dust are lumped together. That would be fine for total AOD but not if you want to get the individual contributions.

Secondly, we are using the ODAER and ODMDUST arrays, which are computed in GeosCore/aerosol_mod.F (or .F90 in 12.8.0). But these arrays hold total aerosol AOD and total dust AOD because they are summed into over all aerosol/dust species. This is where ODAER is summed into over aerosols:

IF ( .not. LUCX .or. ( LUCX .and. (N.LE.NRHAER)) ) THEN
!--------------------------------------------------------
! %%%%%%% Aerosols undergoing hygroscopic growth %%%%%%%
!--------------------------------------------------------
!calculate optics for hyrdophillic aerosol here
!However MDENS in LUT was in g/cm3 not kg/m3 so x1e3
ODAER(I,J,L,IWV,N) = SCALEOD * BXHEIGHT(I,J,L) * 0.75d0 * &
WAERSL(I,J,L,N) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
!Include BC absorption enhancement (xnw, 8/24/15)
IF (N.eq.2) THEN
IF (LBCAE) THEN
BCSCAT_AE = ODAER(I,J,L,IWV,N)*SCALESSA*SSAA(IWV,1,N)
ODAER(I,J,L,IWV,N) = ODAER(I,J,L,IWV,N) * &
( BCAE_1 + SCALESSA*SSAA(IWV,1,N) - &
SCALESSA*SSAA(IWV,1,N)*BCAE_1 )
!now combine with hydrophilic OD as before
BCSCAT_AE = BCSCAT_AE + SSAA(IWV,1,N) * &
0.75d0 * BXHEIGHT(I,J,L) * &
DAERSL(I,J,L,N-1) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
ODAER(I,J,L,IWV,N)= ODAER(I,J,L,IWV,N) + &
(BCAE_2+SSAA(IWV,1,N) - SSAA(IWV,1,N)*BCAE_2) * &
0.75d0 * BXHEIGHT(I,J,L) * &
DAERSL(I,J,L,N-1) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
ELSE
!now combine with hydrophilic OD as before
ODAER(I,J,L,IWV,N)= ODAER(I,J,L,IWV,N) + &
0.75d0 * BXHEIGHT(I,J,L) * &
DAERSL(I,J,L,N-1) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
ENDIF
ENDIF
IF (N.eq.3) THEN
!now combine with hydrophilic OD as before
ODAER(I,J,L,IWV,N)= ODAER(I,J,L,IWV,N) + &
0.75d0 * BXHEIGHT(I,J,L) * &
DAERSL(I,J,L,N-1) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
ENDIF
! Get the AOD contribution from isoprene SOA only (eam, 2014)
IF ( N == 3 .and. Is_ComplexSOA ) THEN
ISOPOD(I,J,L,IWV) = SCALEOD*BXHEIGHT(I,J,L)*0.75d0 &
* ISOAAQ(I,J,L) * QQAA(IWV,1,N) / &
( MSDENS(N) * REAA(1,N) * 1.0D-6 )
ENDIF

and over dust:

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, N )
DO N = 1, NDUST
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
! dust stored in the IDST species bin of LUT variables
ODMDUST(I,J,L,IWV,N) = 0.75e+0_fp * &
State_Met%BXHEIGHT(I,J,L) * &
DUST(I,J,L,N) * QQAA(IWV,N,IDST) / &
( MSDENS(N) * RDAA(N,IDST) * 1.0e-6_fp)
#ifdef RRTMG
!add dust optics to the RT code arrays
!SSA and ASYM copying seems a little redundant...
!will keep this way for uniformity for now but
!possibly could deal with SSA and ASYM in RT module
RTODAER(I,J,L,IWV,NAER+2+N) = ODMDUST(I,J,L,IWV,N)
RTSSAER(I,J,L,IWV,NAER+2+N) = SSAA(IWV,N,IDST)
RTASYMAER(I,J,L,IWV,NAER+2+N) = ASYMAA(IWV,N,IDST)
#endif
ENDDO
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

So long story short, there aren't any arrays that store the individual contributions for aerosol & dust AOD's. You would have to add these and populate them in aerosol_mod and dust_mod.

The planeflight diagnostic is kind of old by now. I am not sure if this ever worked properly, or if it did once and then fell out of step with other model updates since then.

@JFBrewer
Copy link
Author

JFBrewer commented May 8, 2020

Ok, that makes sense, given the values output here. Thank you Bob! I have 2 additional qs:

  1. is the total AOD a column value then? all of these arrays have a L-component - does the AOD change in the course of a single column? Apologies if I'm misunderstanding something here.

  2. As I understand it, bpch outputs are being phased out in 13.0 - does this imply that bpch diagnostics like the Planeflight menu will also be phased out? Asked another way - should I devote more time to working to customize Planeflight output or should I be focusing on learning my way around ObsPack, if that's going to be taking the place of the Planeflight diagnostics?

Thanks,
Jared

@yantosca
Copy link
Contributor

yantosca commented May 8, 2020

The ODAER and ODMDUST are grid box AOD's (since we feed them into FAST-JX). But if you sum them in the vertical you'll get column AOD.

We have no plans to remove Planeflight for the time being, as several people still use it. ObsPack is also a similar diagnostic (and you can use it) but you'd have to convert the input file to the netCDF format that it expects. Also I think ObsPack can save some meteorology and species concs by default but maybe not the AODs.

So your I think your best bet would be to hack Planeflight etc. so that it saves out what you want.

@yantosca
Copy link
Contributor

I am going to close out this issue but feel free to reopen if you have further questions.

@JFBrewer
Copy link
Author

JFBrewer commented May 13, 2020 via email

@msulprizio
Copy link
Contributor

Shixian Zhai (@zsx-GitHub) wrote:

To follow Jared’s finding on the bug of the plane flight AOD diagnostic, I’m attaching my changes (marked with “zhaisx”) on GeosCore/planeflight_mod.F. I made these changes based on GEOS-Chem v12.7.1. To my understanding, the original code incorrectly sums up multiple species and integrates AOD through the entire column for both AODC_tracer and AODB_tracer.

After my changes, planeflight_mod outputs AOD at individual grid and for each tracer in AODC_tracer, and column AOD for each tracer above the aircraft in AODB_tracer. The code can be easily revised to let the model output column AOD in AODC_tracer and column AOD below the aircraft in AODB_tracer as defined on this wiki page.

The reason why I output grid AOD is that grid AOD can be used to derive the extinction coefficient (by dividing by the height of the grid box), which is measured by the campaign I am working on. Also, because the aircraft measured AOD are AOD above the plane in the observation data I use, I output column AOD above the model grid.

Another thing to point out is that in lines 1927-1928 and 2065-2066 in the attached file, ‘ODAER’ were changed to ‘ODMDUST’ in the if statements. This, to my understanding, was a bug in the original code.

Please let me know if I have made any mistakes. Thank you!

@msulprizio msulprizio added this to the 12.9.0 milestone Jun 5, 2020
@msulprizio
Copy link
Contributor

This fix is now merged into our development code in commit 4a3d1b9 and will be included in 12.9.0.

@jennyfisher
Copy link
Contributor

Hi @msulprizio & others,

I've just stumbled across this after encountering the same problem in an older v12 we are running, and I am concerned that the "fix" actually changes the diagnostic outputs so that they no longer match the descriptions given on the wiki (as well as the meaning historically given to these diagnostics).

AODC should be the total column optical depth but is now grid-box AOD, and AODB should be below the aircraft but is now above.

I suggest that we change these back to their original meaning. This would mean:

  • for AODC, adding back in the loops DO LL = 1, State_Grid%NZ and changing L back to LL in the ODAER and OMDUST arrays.
  • for AODB, (I think, less familiar with this one) changing the loop from DO LL = L, State_Grid%NZ to DO LL = 1, L.

If there is a desire for single grid box AOD or above plane AOD it would be easy enough to just add those as additional diagnostics, but the current disconnect is problematic IMO.

Cheers,
Jenny

@zsx-GitHub
Copy link
Contributor

Dear Jenny Fisher & others,

I agree with Jenny that descriptions given on wiki page should match the model to avoid confusion. We can either add additional diagnostics or just edit the wiki. It might be better to use names like AODG and AODA to respectively represent grid-box AOD and above grid AOD.

Best,
Shixian

@jennyfisher
Copy link
Contributor

Hi @msulprizio & @zsx-GitHub --

I discovered another potential issue in the new version: the individual AOD tracers (e.g. AODC_SULF, AODC_BLKC, etc.) do not sum to the total AOD.

This is because when using any simulation that uses UCX (e.g. Standard), there are two stratospheric aerosol tracers. These aren't accounted for explicitly (but were included in the buggy version that looped over all aerosol species).

I think it is likely that users will want to be able to get the total AOD in addition to the individual tracers (as this is what is likely measured). I have a version that adds this as tracer 2007 for total column AODC (haven't messed with AODB as we don't need that for our sims). It's based on a slightly earlier model version, but should be easy to adapt. Let me know if you want it.

-Jenny

@msulprizio
Copy link
Contributor

Hi @jennyfisher. Thanks for following up with the solution. Please do share your fix and we can include it in a future version.

@jennyfisher
Copy link
Contributor

Hi @msulprizio -

Our current version of planeflight_mod.F90 is attached (had to zip for github to accept).

A couple of quick notes:

  1. This version is based on 12.8 code, so before the fixes that described earlier in this thread. We manually implemented the fix to AODC, but in our case this represents total column AOD rather than grid-box specific. This is consistent with the original meaning, implementation, and wiki documentation for AODC. It will be important to make sure that the new AODC_TOT (2007) matches the AODC_XXXX (2001:2006), so either change the others back to total column or change this one to single level (and in that case update the wiki).
  2. We only updated AODC. Haven't touched AODB, so that still has the original bug that spawned this ticket. But it should be easy to start from the updates you already have for that one and just add the total AODB tracer similarly to what I've done for total AODC.
  3. If AODB remains as above-plane AOD as Shixian implemented it, then the wiki page should really be updated.

I'm happy to take a look at an updated version of this code to look for any issues!

-Jenny

planeflight.zip

msulprizio added a commit that referenced this issue Mar 26, 2021
Here we add updates from Jenny Fisher to include total column AOD as
an option for planeflight output. We also add similar updates to include
total AOD below the plane. In addition, we have restored the code so that
it is consistent with the descriptions on the wiki. That is, we loop over
all levels, instead of obtaining AOD for a single level. We also restore
the AODB diagnostics so that they are below plane, not above plane to avoid
confusion.

Original message from Jenny Fisher (Github issue #294):

   I discovered another potential issue in the new version: the individual
   AOD tracers (e.g. AODC_SULF, AODC_BLKC, etc.) do not sum to the total
   AOD.

   This is because when using any simulation that uses UCX (e.g. Standard),
   there are two stratospheric aerosol tracers. These aren't accounted for
   explicitly (but were included in the buggy version that looped over all
   aerosol species).

   I think it is likely that users will want to be able to get the total AOD
   in addition to the individual tracers (as this is what is likely measured).
   I have a version that adds this as tracer 2007 for total column AODC
   (haven't messed with AODB as we don't need that for our sims).

   A couple of quick notes:

   1. This version is based on 12.8 code. We manually implemented the fix to
      AODC, but in our case this represents total column AOD rather than
      grid-box specific. This is consistent with the original meaning,
      implementation, and wiki documentation for AODC. It will be important
      to make sure that the new AODC_TOT (2007) matches the AODC_XXXX
     (2001:2006), so either change the others back to total column or
     change this one to single level (and in that case update the wiki).

   2. We only updated AODC. Haven't touched AODB, so that still has the
      original bug that spawned this ticket. But it should be easy to start
      from the updates you already have for that one and just add the total
      AODB tracer similarly to what I've done for total AODC.

   3. If AODB remains as above-plane AOD as Shixian implemented it, then
      the wiki page should really be updated.

Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>
@msulprizio
Copy link
Contributor

Hi @jennyfisher. We've finally implemented your updates in planeflight_mod.F90 to include total AOD. See commit 6e14981. I also added AODB_TOT following your example. Finally, I reverted AODB_* to represent AOD below the cloud for consistency with what is on the wiki. These updates will be included in 13.1.0.

@jennyfisher
Copy link
Contributor

Great, thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
category: Bug Something isn't working
Projects
None yet
Development

No branches or pull requests

5 participants