geoschem / geos-chem

GEOS-Chem "Science Codebase" repository. Contains GEOS-Chem science routines, run directory generation scripts, and interface code. This repository is used as a submodule within the GCClassic and GCHP wrappers, as well as in other modeling contexts (external ESMs).
http://geos-chem.org
Other
164 stars 156 forks source link

[QUESTION] Curtain along planeflight track #249

Closed yantosca closed 4 years ago

yantosca commented 4 years ago

I am submitting this issue on behalf of Charley Fite (chf18g [at] my.fsu.edu)

I am looking to sample aerosol extinction from a GC simulation along a flight track and I am wanting a vertical curtain along the track. Currently I am using GC version 12.6.0.

I was thinking of using the planeflight module to sample my specific variables at all needed altitudes, for each instance along my flight track, though I’m not sure if this is the recommended way for doing this or if there is a better method?

I don’t see extinction as a variable option in the plane flight wiki, so my method for deriving it is to take the AOD output and divide it by the height thickness of that grid cell. However, I see in the plane flight wiki that the AOD options are for the column integrated value. I wanted to check and make sure whether those are the column integrated values, or if that is the AOD value of that grid cell? For example, if I specified multiple heights at the same lat & lon to get my vertical profile, would all of those have the same AOD value? If so, is there a way to sample the AOD curtain so that I can get my curtain of extinction values?

yantosca commented 4 years ago

Thanks for writing. The planeflight module saves a variable at a single (I,J,L) box corresponding to a flight location, but not a column of points at (I,J). See:

https://github.com/geoschem/geos-chem/blob/c4b68994f9720f246f38a76489bff4913cf489ed/GeosCore/planeflight_mod.F#L1881-L2020

ODAER(I,J,L,IWV,N) ' is the aerosol optical depth at grid-box (I,J,L) for aerosols species N and wavelength IWV. So it is not a column averaged AOD. For the plane flight, this gets averaged into AOD above and AOD below the plane location (seeVARI(V) = VARI(V) + ...`).

VARI(V) is a single (scalar) value at each plane track point. You might be able to modify this into VARI(V,n_lev), where n_lev is the number of vertical grid boxes the column. You would also have to modify the printout to accommodate writing out a column of points. If you have a big plane track, this can make the ASCII file Plane.log pretty big, though. You might even be able to do this better with ObsPack (assuming you can create a netCDF file in the proper format with the coordinate points).

yantosca commented 4 years ago

In case I wasn't clear, you would have to change:

 VARI(V) = VARI(V) + 
 &                ODAER(I,J,LL,IWVSELECT(1,1),ISPC) 

to

 VARI(V,LL) = ODAER(I,J,LL,IWVSELECT(1,1),ISPC)

etc.

chfite commented 4 years ago

Thank you for the response. I figured out something similar to that but am not sure if my thinking is correct. I saw that it loops through all levels with LL and then accumulates those AOD values with VARI(V) = VARI(V) + ...) like you mentioned. Am I correct in thinking that by getting rid of the loop through the levels where it uses LL, and instead I use L to get the value at that specific level along the flight track? Or is there a reason it would still need to loop through all levels with LL?

I made the following edits by adding a new Case:

                  !--------------------------
                  ! Aerosol Optical Depth
                  ! at Particular Grid Cell [unitless]
                  !--------------------------
                  CASE ( 2007:2012 )

                   ! 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
                     ! Get AOD Value at this location and level
                     VARI(V) = ODAER(I,J,L,IWVSELECT(1,1),ISPC)
                    ELSE
                     ! Interpolated using angstrom exponent between
                     ! Closest available wavelengths
                     ! (coefs pre-calculated in CALC_AOD (RD_AOD.F)
                     !catch any zero values before interpolation
                     IF ((ODAER(I,J,L,IWVSELECT(2,1),ISPC).GT.0).AND.
     &                   (ODAER(I,J,L,IWVSELECT(1,1),ISPC).GT.0)) THEN
                          VARI(V) = (ODAER(I,J,L,IWVSELECT(2,1),ISPC)*
     &                    ACOEF_WV(1)**(BCOEF_WV(1)*
     &                    LOG(ODAER(I,J,L,IWVSELECT(1,1),ISPC)/
     &                    ODAER(I,J,L,IWVSELECT(2,1),ISPC))))
                     ENDIF
                    ENDIF
                   ENDDO

                   !now add in the dust
                   IF ( .not. LINTERP ) THEN
                    DO ISPC = 1, NDUST
                     ! Get AOD Value at this location and level
                     VARI(V) = ODMDUST(I,J,L,IWVSELECT(1,1),ISPC)
                    ENDDO
                   ELSE
                    ! Interpolated using angstrom exponent between
                    ! Closest available wavelengths
                    ! (coefs pre-calculated in CALC_AOD (RD_AOD.F)
                    !catch any zero values before interpolation
                    DO ISPC = 1, NDUST
                     IF ((ODAER(I,J,L,IWVSELECT(2,1),ISPC).GT.0).AND.
     &                   (ODAER(I,J,L,IWVSELECT(1,1),ISPC).GT.0)) THEN
                          VARI(V) = (ODMDUST(I,J,L,IWVSELECT(2,1),ISPC)*
     &                    ACOEF_WV(1)**(BCOEF_WV(1)*
     &                    LOG(ODMDUST(I,J,L,IWVSELECT(1,1),ISPC)/
     &                    ODMDUST(I,J,L,IWVSELECT(2,1),ISPC))))
                     ENDIF
                    ENDDO
                   ENDIF
yantosca commented 4 years ago

As is my understanding (I,J,L) is the grid box where the plane is. If you want to go above or below that you need the loop over LL.

chfite commented 4 years ago

In order to get my vertical curtain above and below the plane, I am specifying all of those heights in the planeflight input file, giving a point for each height above and below the plane along the flight track. So, I think for each individual point/height that I feed in along the track I would want to sample at (I,J,L) with L being a height that could be above, below, or where the plane is. That was my reasoning for not looping through levels with LL.

yantosca commented 4 years ago

Ah, OK, now I understand. Try that and let us know if it works.

chfite commented 4 years ago

My methodology for getting the AOD in Planeflight at a particular grid-cell (as opposed to column integrated) worked, as I was able to get an output of AOD above and below my plane in a curtain and it looked as it should.

One new question/issue I had though was it appears that no matter what AOD type you choose in planeflight_mod, it will sample the same AOD values. I think this is because it loops through NAER with ISPC and gets an AOD value that corresponds to ISPC. Does ISPC represent which type of AOD species it is in ODAER? If so, that would mean that by including the loop with ISPC in it is getting the total AOD by combing all species. Like below.

! Loop over RH bins DO ISPC= 1, NAER !Get AOD Value at this location and level VARI(V) = VARI(V) +ODAER(I,J,L,IWVSELECT(1,1),ISPC) Keep in mind in the above code example I removed the loop through LL because I want AOD at a specific level L, not the column integrated.

I think if one was wanting to get the AOD at a particular level AND for just a particular species (not just total AOD), then the code would have to look something below by removing the loop of ISPC and specify the particular ISPC value that corresponds to the AOD species you'd want. You'd also have to remove VARI(V) + ... portion within the loop. VARI(V) = ODAER(I,J,L,IWVSELECT(1,1),ISPC) In this case how would I know which value for ISPC would represent the specific AOD species I want (e.g., BC, OC, etc.)?

Or perhaps my understanding of ISPC is here is incorrect? If so, how how is planeflight set up to sample a specific AOD species of choice?

yantosca commented 4 years ago

Thanks for writing. If you look in aerosol_mod.F90 you'll see this code:

https://github.com/geoschem/geos-chem/blob/c4b68994f9720f246f38a76489bff4913cf489ed/GeosCore/aerosol_mod.F#L1412-L1444

which gives the mapping of the slots of ODAER. 1 = SO4 2 = BCPI 3= OCPI (or POA if complex SOA) 4 = SALA 5 = SALC

For UCX-based simulations there are 2 more slots:

6 = SSA/STS 7 = NAT/ice PSC

but for tropchem there are only 5 slots.

I agree -- this isn't really well documented. You have to dig into the code to find this. Aerosol_mod is a bit of a mess. I think in the future we might try to clean it up a bit but that might be a ways off.

Hope this helps!

yantosca commented 4 years ago

I am going to close out this issue. Feel free to reopen.