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

Tagged o3 simulation- tagged stratospheric ozone producing values of 0 throughout the troposphere #2028

Open ncolombi opened 10 months ago

ncolombi commented 10 months ago

Name and Institution (Required)

Name: Nadia Colombi Institution: Harvard

Confirm you have reviewed the following documentation

Description of your issue or question

I ran the GC version 14.2.0 tagged ozone simulation, and it finished with no errors. I allowed 5 months of spinup for the run, and the free tropospheric ozone (O3_VV) looks normal (figure 1 below), however, the O3strat_VV diagnostic is showing values of 0 for everywhere in the atmosphere apart from the upper stratosphere (figures 2, 3, 4). I am using ProdLoss files from a previous run so that regional emissions can be included. I ran a separate simulation out of the box with 1.5 years of spinup to make sure the ProdLoss files weren't the issue, but ran into the same problem. My HEMCO_Config.rc and geoschem_config.yml files are attached.

Figure 1

Free tropospheric (800-400 hPa) ozone (O3_VV in ppb) averaged for May-June Screen Shot 2023-10-30 at 2 31 50 PM

Figure 2

Tagged stratosphere (O3strat_VV) ozone in ppb for the free troposphere for May-June Screen Shot 2023-10-30 at 2 30 33 PM (1)

Figure 3

Tagged stratospheric (O3strat_VV) ozone profile (pressure vs. ozone in ppb) Screen Shot 2023-11-01 at 2 15 25 PM

Figure 4

Zonal tagged stratosphere ozone (O3strat_VV) plot for May-June

Screen Shot 2023-11-01 at 2 06 31 PM

Please provide as much detail as possible. Always include the GEOS-Chem version number and any relevant configuration and log files.

files.zip

yantosca commented 10 months ago

Hi @ncolombi. I took your geoschem_config.yml and HEMCO_Config.rc files and did an "out-of-the-box" tagged O3 simulation using your generated prod/loss rates. I came to the same conclusion as you:

nadia1

nadia2

I think I see what is happening here. In the module GeosCore/tagged_o3_mod.F90 we have:

    ! Loop over the # of advected species
    !$OMP PARALLEL DO                                                        &
    !$OMP DEFAULT( SHARED                                                   )&
    !$OMP PRIVATE( I, J, L, N, LO3_kg, LO3_kgps, PO3_kg, PO3_kgps           )&
    !$OMP COLLAPSE( 3                                                       )&
    !$OMP SCHEDULE( DYNAMIC, 4                                              )
    DO L = 1, State_Grid%NZ
    DO J = 1, State_Grid%NY
    DO I = 1, State_Grid%NX

       !=====================================================================
       ! Convert P(O3) and L(O3) from [molec/cm3/s] -> [kg]
       ! cf. https://github.com/geoschem/geos-chem/issues/1109 by Xingpei Ye
       !=====================================================================

       ! P(O3) in [kg/s] and [kg]
       PO3_kgps = PO3_hco(I,J,L)              & ! in molec/cm3/s
                * molcm3_to_kgm3              & ! molec/cm3/s -> kg/m3/s
                * State_Met%AIRVOL(I,J,L)       ! kg/m3/s     -> kg/s
       PO3_kg   = PO3_kgps * dtChem             ! kg/s        -> kg

       ! L(O3) in [kg/s] and [kg]
       LO3_kgps = LO3_hco(I,J,L)              & ! in molec/cm3/s
                * molcm3_to_kgm3              & ! molec/cm3/s -> kg/m3/s
                * State_Met%AIRVOL(I,J,L)       ! kg/m3/s     -> kg/s
       LO3_kg  = LO3_kgps * dtchem              ! kg/s        -> kg

       ! Prevent denormal values
       PO3_kg   = MAX( PO3_kg,   1.0e-30_fp )
       PO3_kgps = MAX( PO3_kgps, 1.0e-30_fp )
       LO3_kg   = MAX( LO3_kg,   1.0e-30_fp )
       LO3_kgps = MAX( LO3_kgps, 1.0e-30_fp )

       ... etc skipped ...

       !=====================================================================
       ! Apply chemical production of ozone (only where it is produced)
       !=====================================================================

       ! Add P(O3) [kg] to the total O3 species
       Spc(id_O3)%Conc(I,J,L) = Spc(id_O3)%Conc(I,J,L) + PO3_kg

       ! Add P(O3) [kg] to the stratospheric O3 species
       IF ( State_Met%InStratosphere(I,J,L) ) THEN
          Spc(id_O3Strat)%Conc(I,J,L) = Spc(id_O3Strat)%Conc(I,J,L) + PO3_kg
       ENDIF

I think the State_Met%InStratosphere(I,J,L) is excluding the mesosphere... this is only TRUE if you are in the stratosphere but not the mesosphere. Try replacing this with State_Met%InStratMeso(I,J,L), which should be TRUE if you are anywhere above the tropopause.

yantosca commented 10 months ago

Also tagging @msulprizio

yantosca commented 10 months ago

@ncolombi: We can add a fix for this in the dev/no-diff-to-benchmark branch, which will be mergeable at any time.