About photolysis_mod.F90 #2641

ChenBHXMU · 2024-12-16T08:56:49Z

Your name

Baihua Chen

Your affiliation

Xiamen University

Please provide a clear and concise description of your question or discussion topic.

I hope this message finds you well!

Currently, sea salt only includes SALC and SALA. I am researching how to add two mode for sea salt. However, I encountered the following issues while optimizing the code. I have tried my best to resolve them, but I still don't know how to modify the photolysis_mod.F90 file. I don't understand the specific values in the code IND = (/22,29,36,43,50/). How were 22, 29, 36, 43, 50 obtained? How can I optimize this code at the moment?

Thank you for your help!

Best regards,
Baihua Chen

lizziel · 2024-12-16T19:56:21Z

Hi @ChenBHXMU, the indexes there have been hard-coded since 2013, with a comment that they are from aerosol_mod.F (now aerosol_mod.F90). That reference is due to an expected order of the hygroscopic growth aerosols that was defined in that module, specifically here:

geos-chem/GeosCore/aerosol_mod.F90

Lines 2481 to 2497 in 6bbd81f

    
           ! Set the mapping to the ordering of aerosol densities in RD_AOD 
        
           SELECT CASE ( TRIM(SpcInfo%Name) ) 
        
           CASE ( 'SO4' ) 
        
              Map_NRHAER(N) = 1 
        
           CASE ( 'BCPI' ) 
        
              Map_NRHAER(N) = 2 
        
           CASE ( 'OCPI', 'POA1' ) 
        
              Map_NRHAER(N) = 3 
        
           CASE ( 'SALA' ) 
        
              Map_NRHAER(N) = 4 
        
           CASE ( 'SALC' ) 
        
              Map_NRHAER(N) = 5 
        
           CASE DEFAULT 
        
              ErrMsg = 'WARNING: aerosol diagnostics not defined' // & 
        
                       ' for NRHAER greater than 5!' 
        
              CALL GC_ERROR( ErrMsg, RC, 'Init_Aerosol in aerosol_mod.F90' ) 
        
           END SELECT

NRHAER is referring to the number of the hygroscopic growth aerosols in the model defined in CMN_SIZE_mod.F90.

geos-chem/Headers/CMN_SIZE_mod.F90

Lines 45 to 46 in 6bbd81f

    
           ! Number of aerosols undergoing hygroscopic growth 
        
           INTEGER, PARAMETER :: NRHAER = 5

So where does IND = (/22,29,36,43,50/) come from? The key is to understand that IND is used to set elements in array MIEDX for hygroscopic growth aerosols, and that MIEDX contains indexes used to access optical properties from file FJX_spec-aer.dat. MIEDX contains mappings to the file for other things as well. For example, the first values are set for clouds and dust here.

geos-chem/GeosCore/photolysis_mod.F90

Lines 850 to 862 in 6bbd81f

    
           ! Clouds 
        
           MIEDX(1)  =  3   !  Black carbon absorber 
        
           MIEDX(2)  = 10   !  Water Cloud (Deirmenjian 8 micron) 
        
           MIEDX(3)  = 14   !  Irregular Ice Cloud (Mishchenko) 
        
           ! Dust 
        
           MIEDX(4)  = 15   !  Mineral Dust  .15 micron    (rvm, 9/30/00) 
        
           MIEDX(5)  = 16   !  Mineral Dust  .25 micron    (rvm, 9/30/00) 
        
           MIEDX(6)  = 17   !  Mineral Dust  .4  micron    (rvm, 9/30/00) 
        
           MIEDX(7)  = 18   !  Mineral Dust  .8  micron    (rvm, 9/30/00) 
        
           MIEDX(8)  = 19   !  Mineral Dust 1.5  micron    (rvm, 9/30/00) 
        
           MIEDX(9)  = 20   !  Mineral Dust 2.5  micron    (rvm, 9/30/00) 
        
           MIEDX(10) = 21   !  Mineral Dust 4.0  micron    (rvm, 9/30/00)

You can take a look at your Cloud-J dat-file FJX_spec-aer.dat and look for index 03. It will look like this, with a comment that it corresponds to soot. You can try to find the others listed above and see how they map to the file.

  03|       ABSRB|  0.039  1.000| fully absorbing 'soot', wavelength indep.
 200   1.0000 0.0000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 300   1.0000 0.0000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 400   1.0000 0.0000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 600   1.0000 0.0000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 999   1.0000 0.0000  0.000  0.000  0.000  0.000  0.000  0.000  0.000

In total, the MIEDX array has 3 entries for clouds, 7 entries for dust, 7 entries for each of the 5 hygroscopic growth aerosols used in the model, and 2 entries for stratospheric aerosols. This is because we have optical properties for each of those stored in FJX_spec-aer.dat. The IND array stores the first index in the dat-file for each of the hygroscopic growth aerosols. The values are in increments of 7 because the dat-file stores properties for seven different relative humidity ranges per aerosol. If you look for index 22 in the file you will find 22| S00(dar)| 0.121 1.700| Trop sulfate, RH=00. The is the first relative humidity bin for tropospheric sulfate. Notice that order of the aerosols in the dat-file is the same as the order of the aerosols in aerosol_mod.F90 when Map_NRHAER is defined (see earlier code snippet).

The mappings in MIEDX for hygroscopic growth aerosols are set like this:

geos-chem/GeosCore/photolysis_mod.F90

Lines 864 to 869 in 6bbd81f

    
           ! Aerosols 
        
           DO I=1,NRHAER 
        
              DO J=1,NRH 
        
                 MIEDX(10+((I-1)*NRH)+J)=IND(I)+J-1 
        
              ENDDO 
        
           ENDDO

Here is where it is important to understand where these constants are coming from. All are stored in file CMN_SIZE_mod.F90, and these values are important if you are adding aerosols in the model.

geos-chem/Headers/CMN_SIZE_mod.F90

Lines 42 to 55 in 6bbd81f

    
           ! Number of FAST-J aerosol size bins (rvm, bmy, 11/15/01) 
        
           INTEGER, PARAMETER :: NDUST = 7 
        
           ! Number of aerosols undergoing hygroscopic growth 
        
           INTEGER, PARAMETER :: NRHAER = 5 
        
           ! Number of stratospheric aerosols (SDE 04/17/13) 
        
           INTEGER, PARAMETER :: NSTRATAER = 2 
        
           ! Number of other aerosol categories, include stratospheric aerosols 
        
           INTEGER, PARAMETER :: NAER = NRHAER + NSTRATAER 
        
           ! NRH -- number of relative humidity bins (rvm, bmy, 2/27/02) 
        
           INTEGER,  PARAMETER :: NRH = 5

You should only need to expand IND if you want to account for extinction due to the new aerosols in photolysis. If you do, then you will also need to add entries to FJX_spec-aer.dat and update setting MIEDX accordingly.

I hope this helps!

ChenBHXMU · 2024-12-17T09:03:42Z

Dear @lizziel ,

Thank you for your patient guidance. However, the FJX_spec-aer.dat file seem be no longer accessible. I have found the FJX_spec.dat file, but its content does not match that of the FJX_spec-aer.dat. I cann't locate the FJX_spec-aer.dat file. What steps should I take to resolve this issue?

Thank you so much!

Best regards,
Baihua Chen

lizziel · 2024-12-17T15:39:37Z

Hi @ChenBHXMU, the tables directory in Cloud-J is used for Cloud-J standalone and not GEOS-Chem. The directory containing the file is configured within geoschem_config.yml.

ChenBHXMU · 2024-12-18T08:51:14Z

Dear @lizziel ,

I'm glad to receive your message. However, I still can't find the FJX_spec-aer.dat file after searching my computer globally. A friend pointed out that the file might be ExtData/CHEM_INPUTS/FAST_JX/v2019-06/jv_spec_mie.dat, but I also couldn't find this file. I think it might be due to version differences. My version is 14.4.3.

Could you provide further assistance? Thank you for your help.

Best regards,
Baihua Chen

lizziel · 2024-12-18T14:52:38Z

Ah, I think the confusion is the filename is actually FJX_scat-aer.dat.

lizziel · 2024-12-18T14:58:35Z

Please note you can search for strings in the code using grep as a way to find things. For example, if you do this at the highest level of GEOS-Chem Classic 14.4.3 code:

grep -r "aer.dat" *

Then you will get this, which shows the correct name of the file:

src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:    ! FJX_spec.dat (RD_XXX), FJX_scat-aer.dat (RD_MIE), and 
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:    ! Relative humidities in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                ! relative humidity entries in FJX_spec-aer.dat when computing optical
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                ! one for each relative humidity entry in FJX_spec-aer.dat. Values are
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:                   K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat
src/GEOS-Chem/GeosCore/cldj_interface_mod.F90:    filePath = TRIM( Input_Opt%CloudJ_Dir ) // 'FJX_scat-aer.dat'
src/Cloud-J/tables/FJX_scat-aer.dat:(FJX_scat-aer.dat) Aerosol scatter phase fns (v 7.3+, new format, old data)
src/Cloud-J/src/Core/cldj_init_mod.F90:      call RD_MIE(AMIROOT,JXUNIT,TRIM(DATADIR)//'FJX_scat-aer.dat',rc)
src/Cloud-J/src/Core/cldj_cmn_mod.F90:! Variables in file 'FJX_scat-aer.dat' (RD_MIE)
src/Cloud-J/CHANGELOG.md:- Changed FJX_scat-aer.dat format to expect 1 character greater width for R-eff and rho
src/Cloud-J/CHANGELOG.md:- Changed FJX_scat-aer.dat format to expect 2 characters greater width for QAA (data column 2)
src/Cloud-J/CHANGELOG.md:- Changed FJX_scat-aer.dat format to expect 1 character greater width for PAA (data columns 3-9)
src/Cloud-J/docs/CloudJ77_notes.txt:    6,760 FJX_scat-aer.dat
src/Cloud-J/docs/CloudJ77_notes.txt:       6,760 FJX_scat-aer.dat

Getting familiar with this method of searching will be essential if you plan on doing a major code update to the model.

ChenBHXMU · 2024-12-20T08:12:03Z

Dear @lizziel ,

Thank you for your patient guidance. However, I still have some issues that I cannot resolve. I would like to further subdivide the Fine mode and the Coarse mode. Split the Fine mode into two, and also split the Coarse mode into two. The main modifications are as follows:

geos-chem/Headers/CMN_SIZE_mod.F90：
! Number of aerosols undergoing hygroscopic growth Initial NRHAER = 5; by cbh INTEGER, PARAMETER :: NRHAER = 7

GEOS-Chem/GeosCore/photolysis_mod.F90:
IND = (/22,29,36,43,43,50,50/)

GEOS_Chem/GeosCore/aerosol_mod.F90:
` DO N = 1, NRHAER
! Get the species database index from the species database
! mapping array for hygroscopic growth species
SpcID = State_Chm%Map_HygGrth(N)

   ! Point to the Species Database entry for species N
   SpcInfo => State_Chm%SpcData(SpcID)%Info

   ! Set the mapping to the ordering of aerosol densities in RD_AOD
   SELECT CASE ( TRIM(SpcInfo%Name) )
   CASE ( 'SO4' )
      Map_NRHAER(N) = 1
   CASE ( 'BCPI' )
      Map_NRHAER(N) = 2
   CASE ( 'OCPI', 'POA1' )
      Map_NRHAER(N) = 3
   CASE ( 'SALA')
      Map_NRHAER(N) = 4
   case ( 'SALA2' )
      Map_NRHAER(N) = 5
   CASE ( 'SALC' )
      Map_NRHAER(N) = 6
   CASE ( 'SALC2' )
      Map_NRHAER(N) = 7
    
   CASE DEFAULT
      ErrMsg = 'WARNING: aerosol diagnostics not defined' // &
               ' for NRHAER greater than 5!'
      CALL GC_ERROR( ErrMsg, RC, 'Init_Aerosol in aerosol_mod.F90' )
   END SELECT

   ! Free pointer
   SpcInfo => NULL()

ENDDO`

GEOS-Chem/GeosCore/cldj_interface_mod.F90:
` !----------------------------------------------------
! Seasalt [dry molec/cm3] -> [wet g/m2]
!----------------------------------------------------

            IF ( Input_Opt%LSSALT ) THEN

               !----------------------------------------------------
               ! Accumulation mode seasalt
               !----------------------------------------------------
               !WRITE(*,*) 'myNRHAER:', NRHAER*(SALA_ind-1)
               ! Get indexes to optical property LUT
               S_rh0 = 3 + NDUST + NRHAER*(SALA_ind-1) + 1  ! SALA index for RH=0 in NDXAER
               S_rhx = S_rh0 + RH_ind - 1  ! SALA index for this RH
               K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
               K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat

               ! Get interpolated effective radius and extinction for RH in this grid box
               IF ( RH_ind == NRH ) THEN
                  RAA_eff = RAA(K_rhx)
                  QAA_eff = QAA(ind_1000,K_rhx)
               ELSE
                  FRAC = ( State_Met%RH(I,J,L) - RH_lut(RH_ind) ) &
                       / ( RH_lut(RH_ind+1) - RH_lut(RH_ind) )
                  RAA_eff = RAA(K_rhx) + FRAC * ( RAA(K_rhx+1) - RAA(K_rhx) )
                  QAA_eff = QAA(ind_1000,K_rhx) + FRAC * ( QAA(ind_1000,K_rhx+1) - QAA(ind_1000,K_rhx) )
               ENDIF
               dry_to_wet_factor = ( RAA_eff / RAA(K_rh0) )**3
               R_interp_factor = RAA_eff / RAA(K_rhx)
               Q_interp_factor = QAA_eff / QAA(ind_1000,K_rhx)

               ! Set concentration, converting [dry kg/m3] -> [wet g/m2]
               AERSP(L,S_rhx) = State_Chm%AerMass%SALA(I,J,L) * BoxHt * 1.d3 * dry_to_wet_factor &
                    * Q_interp_factor / R_interp_factor
                    
               
               !----------------------------------------------------
               ! Accumulation mode seasalt-2
               !----------------------------------------------------
               !WRITE(*,*) 'myNRHAER:', NRHAER*(SALA_ind2-1)
               ! Get indexes to optical property LUT
               S_rh0 = 3 + NDUST + NRHAER*(SALA_ind2-1) + 1  ! SALA2 index for RH=0 in NDXAER
               S_rhx = S_rh0 + RH_ind - 1  ! SALA index for this RH
               K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
               K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat

               ! Get interpolated effective radius and extinction for RH in this grid box
               IF ( RH_ind == NRH ) THEN
                  RAA_eff = RAA(K_rhx)
                  QAA_eff = QAA(ind_1000,K_rhx)
               ELSE
                  FRAC = ( State_Met%RH(I,J,L) - RH_lut(RH_ind) ) &
                       / ( RH_lut(RH_ind+1) - RH_lut(RH_ind) )
                  RAA_eff = RAA(K_rhx) + FRAC * ( RAA(K_rhx+1) - RAA(K_rhx) )
                  QAA_eff = QAA(ind_1000,K_rhx) + FRAC * ( QAA(ind_1000,K_rhx+1) - QAA(ind_1000,K_rhx) )
               ENDIF
               WRITE(*,*) "myRAA(K_rh0):", RAA(K_rh0)
               dry_to_wet_factor = ( RAA_eff / RAA(K_rh0) )**3
               R_interp_factor = RAA_eff / RAA(K_rhx)
               Q_interp_factor = QAA_eff / QAA(ind_1000,K_rhx)

               ! Set concentration, converting [dry kg/m3] -> [wet g/m2]
               AERSP(L,S_rhx) = State_Chm%AerMass%SALA2(I,J,L) * BoxHt * 1.d3 * dry_to_wet_factor &
                    * Q_interp_factor / R_interp_factor
               
               

               !----------------------------------------------------
               ! Coarse seasalt
               !----------------------------------------------------

               ! Get indexes to optical property LUT
               S_rh0 = 3 + NDUST + NRHAER*(SALC_ind-1) + 1  ! SALC index for RH=0 in NDXAER
               S_rhx = S_rh0 + RH_ind - 1  ! SALC index for this RH
               K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
               K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat

               ! Get interpolated effective radius and extinction for RH in this grid box
               IF ( RH_ind == NRH ) THEN
                  RAA_eff = RAA(K_rhx)
                  QAA_eff = QAA(ind_1000,K_rhx)
               ELSE
                  FRAC = ( State_Met%RH(I,J,L) - RH_lut(RH_ind) ) &
                       / ( RH_lut(RH_ind+1) - RH_lut(RH_ind) )
                  RAA_eff = RAA(K_rhx) + FRAC * ( RAA(K_rhx+1) - RAA(K_rhx) )
                  QAA_eff = QAA(ind_1000,K_rhx) + FRAC * ( QAA(ind_1000,K_rhx+1) - QAA(ind_1000,K_rhx) )
               ENDIF
               dry_to_wet_factor = ( RAA_eff / RAA(K_rh0) )**3
               R_interp_factor = RAA_eff / RAA(K_rhx)
               Q_interp_factor = QAA_eff / QAA(ind_1000,K_rhx)

               ! Set concentration, converting [dry molec/cm3] -> [wet g/m2]
               AERSP(L,S_rhx) = State_Chm%AerMass%SALC(I,J,L) * BoxHt  * 1.d3 * dry_to_wet_factor &
                    * Q_interp_factor / R_interp_factor
                    
               
               !----------------------------------------------------
               ! Coarse seasalt-2
               !----------------------------------------------------

               ! Get indexes to optical property LUT
               S_rh0 = 3 + NDUST + NRHAER*(SALC_ind2-1) + 1  ! SALC index for RH=0 in NDXAER
               S_rhx = S_rh0 + RH_ind - 1  ! SALC index for this RH
               K_rh0 = NDXAER(L,S_rh0)     ! index for RH=0 in FJX_spec-aer.dat
               K_rhx = NDXAER(L,S_rhx)     ! index for this RH in FJX_spec-aer.dat

               ! Get interpolated effective radius and extinction for RH in this grid box
               IF ( RH_ind == NRH ) THEN
                  RAA_eff = RAA(K_rhx)
                  QAA_eff = QAA(ind_1000,K_rhx)
               ELSE
                  FRAC = ( State_Met%RH(I,J,L) - RH_lut(RH_ind) ) &
                       / ( RH_lut(RH_ind+1) - RH_lut(RH_ind) )
                  RAA_eff = RAA(K_rhx) + FRAC * ( RAA(K_rhx+1) - RAA(K_rhx) )
                  QAA_eff = QAA(ind_1000,K_rhx) + FRAC * ( QAA(ind_1000,K_rhx+1) - QAA(ind_1000,K_rhx) )
               ENDIF
               dry_to_wet_factor = ( RAA_eff / RAA(K_rh0) )**3
               R_interp_factor = RAA_eff / RAA(K_rhx)
               Q_interp_factor = QAA_eff / QAA(ind_1000,K_rhx)

               ! Set concentration, converting [dry molec/cm3] -> [wet g/m2]
               AERSP(L,S_rhx) = State_Chm%AerMass%SALC2(I,J,L) * BoxHt  * 1.d3 * dry_to_wet_factor &
                    * Q_interp_factor / R_interp_factor

            ENDIF

         ENDIF`

` ! Aerosol indexes (must match mapping set in RD_AOD)

SO4_ind  = 1
BC_ind   = 2
OC_ind   = 3
SALA_ind = 4
SALA_ind2 = 5
SALC_ind = 6
SALC_ind2 = 7`

The issue arises from: R_interp_factor = RAA_eff / RAA(K_rhx), where RAA(K_rhx) could be 0.

Thank you for your help, looking forward to your reply.

Best regards,
Baihua Chen

lizziel · 2024-12-20T17:27:30Z

Hi Baihua, if you are reusing the same optical properties for the new species, e.g. same for both bins of coarse, I would think you could add the concentrations together to use for photolysis. This would avoid the hassle of expanding arrays used in photolysis based on the new species. If this is not a viable option for your research then I suggest that you add prints to figure out why there is no value for RAA for certain values of K_rhx. Given you are changing a core size in the model I also suggest you very carefully review all arrays that use NRHAER and any other indexes derived from it to make sure they are correct. Since this is not a research project the GEOS-Chem Support Team is involved in we can only provide limited coding help.

ChenBHXMU added the category: Question Further information is requested label Dec 16, 2024

lizziel self-assigned this Dec 16, 2024

yantosca added topic: Aerosols Related to aerosol species in GEOS-Chem topic: Photolysis Related to photolyis rate computations labels Dec 17, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

About photolysis_mod.F90 #2641

About photolysis_mod.F90 #2641

ChenBHXMU commented Dec 16, 2024

lizziel commented Dec 16, 2024

ChenBHXMU commented Dec 17, 2024

lizziel commented Dec 17, 2024

ChenBHXMU commented Dec 18, 2024

lizziel commented Dec 18, 2024

lizziel commented Dec 18, 2024

ChenBHXMU commented Dec 20, 2024 •

edited

Loading

lizziel commented Dec 20, 2024

About photolysis_mod.F90 #2641

About photolysis_mod.F90 #2641

Comments

ChenBHXMU commented Dec 16, 2024

Your name

Your affiliation

Please provide a clear and concise description of your question or discussion topic.

lizziel commented Dec 16, 2024

ChenBHXMU commented Dec 17, 2024

lizziel commented Dec 17, 2024

ChenBHXMU commented Dec 18, 2024

lizziel commented Dec 18, 2024

lizziel commented Dec 18, 2024

ChenBHXMU commented Dec 20, 2024 • edited Loading

lizziel commented Dec 20, 2024

ChenBHXMU commented Dec 20, 2024 •

edited

Loading