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

Add fixes for transport tracers simulation #1890

Merged
merged 5 commits into from
Jul 28, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
- Added GCHP run-time option in GCHP.rc to choose dry or total pressure in GCHP advection
- Added GCHP run-time option in GCHP.rc to correct native mass fluxes for humidity
- Added new tracer_mod.F90 containing subroutines for applying sources and sinks for the TransportTracer simulation
- Added new species to the TransportTracer simulation: aoa (replaces CLOCK), aoa_bl, aoa_nh, st80_25, stOX
- Added new species to the TransportTracer simulation: aoa (replaces CLOCK), aoa_bl, aoa_nh, st80_25
- Added GEOS-IT and GEOSIT as allowable meteorology source options in geoschem_config.yml

### Changed
Expand Down
193 changes: 50 additions & 143 deletions GeosCore/tracer_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
REAL(fp) :: Local_Tally
REAL(fp) :: Total_Area
REAL(fp) :: Total_Spc
REAL(fp) :: Flux

! SAVEd scalars
LOGICAL, SAVE :: First = .TRUE.
Expand All @@ -113,12 +114,8 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
CHARACTER(LEN=63) :: OrigUnit

! Arrays
REAL(fp) :: Flux(State_Grid%NX,State_Grid%NY)
REAL(fp) :: Mask(State_Grid%NX,State_Grid%NY,State_Grid%NZ)

! Pointers to fields in the HEMCO data structure
REAL(f4), POINTER :: O3(:,:,:) => NULL()

! Objects
TYPE(Species), POINTER :: SpcInfo

Expand All @@ -140,7 +137,6 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
Local_Tally = 0.0_fp
Total_Area = 0.0_fp
Total_Spc = 0.0_fp
Mask = 1.0_fp
Flux = 0.0_fp

#if defined( MODEL_GEOS ) || defined( MODEL_GCHP )
Expand All @@ -154,7 +150,7 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, &
'v/v dry', RC, OrigUnit=OrigUnit )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Unit conversion error (kg/kg dry -> v/v dry)'
ErrMsg = 'Unit conversion error (kg -> v/v dry)'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
Expand All @@ -172,6 +168,9 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
! Skip Rn-Pb-Be tracers for now
IF ( SpcInfo%Is_RadioNuclide ) CYCLE

! Initialize mask for each species
Mask = 1.0_fp

! Skip if this species does not have a source or if the
! source is handled by HEMCO
IF ( TRIM(SpcInfo%Src_Mode) == 'none' .or. &
Expand All @@ -180,7 +179,7 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
IF ( First ) THEN

!------------------------------------------------------------------
! Convert Src_Value to v/v or to seconds for aoa tracers
! Convert Src_Value to v/v
!------------------------------------------------------------------
IF(TRIM(SpcInfo%Src_Units) == 'ppmv' ) THEN

Expand All @@ -194,6 +193,9 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &

SpcInfo%Src_Value = SpcInfo%Src_Value * 1.0E-12

!------------------------------------------------------------------
! Convert timesteps (s) to expected units (e.g. days)
!------------------------------------------------------------------
ELSE IF ( TRIM(SpcInfo%Src_Units) == 'timestep' ) THEN

IF ( TRIM(SpcInfo%Units) == 'hours' ) THEN
Expand Down Expand Up @@ -221,33 +223,6 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &

ENDIF

!---------------------------------------------------------------------
! Get model fields (if needed)
!---------------------------------------------------------------------
IF ( TRIM(SpcInfo%Src_Mode) == 'model_field' ) THEN

! For stOX tracer get O3 [v/v] from HEMCO
IF ( TRIM(SpcInfo%Name ) == 'stOX' ) THEN
CALL HCO_GetPtr( HcoState, 'GLOBAL_O3', O3, RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Cannot get pointer to GLOBAL_O3!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF

ELSE

ErrMsg = 'Src_Mode: model_field not currently supported for' // &
TRIM( SpcInfo%Name ) // &
'Please modify species_database.yml or add this ' // &
'capability in tracer_mod.F90.'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN

ENDIF

ENDIF

!---------------------------------------------------------------------
! Determine mask
!---------------------------------------------------------------------
Expand Down Expand Up @@ -350,18 +325,16 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX

IF ( SpcInfo%Src_Add ) THEN

! Add source
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) + &
( SpcInfo%Src_Value * Mask(I,J,L) )

ELSE
IF ( Mask(I,J,L) > 0 ) THEN

! Replace value
State_Chm%Species(N)%Conc(I,J,L) = &
( SpcInfo%Src_Value * Mask(I,J,L) )
IF ( SpcInfo%Src_Add ) THEN
! Add source
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) + SpcInfo%Src_Value
ELSE
! Replace value
State_Chm%Species(N)%Conc(I,J,L) = SpcInfo%Src_Value
ENDIF

ENDIF

Expand All @@ -370,37 +343,15 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
ENDDO
!$OMP END PARALLEL DO

ELSE IF ( TRIM(SpcInfo%Src_Mode) == 'model_field' ) THEN

! Special handling for stOX tracer
IF ( TRIM(SpcInfo%Name) == 'stOX' ) THEN

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED )&
!$OMP PRIVATE( I, J, L )&
!$OMP COLLAPSE( 3 )
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
! Replace value
State_Chm%Species(N)%Conc(I,J,L) = O3(I,J,L) * Mask(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF

ELSE IF ( TRIM(SpcInfo%Src_Mode) == 'maintain_mixing_ratio' ) THEN

! To distrubute tracer uniformly on the surface, compute the
! total area
! total area [m2]
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX

! Accumulate the total area [m2]
Local_Tally = Local_Tally + &
( State_Grid%Area_M2(I,J) * Mask(I,J,1) )

IF ( Mask(I,J,1) > 0 ) THEN
Local_Tally = Local_Tally + State_Grid%Area_M2(I,J)
ENDIF
ENDDO
ENDDO

Expand All @@ -418,10 +369,12 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
Local_Tally = Local_Tally &
IF ( Mask(I,J,L) > 0 ) THEN
Local_Tally = Local_Tally &
+ ( SpcInfo%Src_Value - State_Chm%Species(N)%Conc(I,J,L) ) &
* ( State_Met%AIRNUMDEN(I,J,L) / AVO ) &
* State_Met%AIRVOL(I,J,L)
ENDIF
ENDDO
ENDDO
ENDDO
Expand All @@ -434,13 +387,13 @@ SUBROUTINE Tracer_Source_Phase( Input_Opt, State_Chm, State_Grid, &
#endif

! Compute flux [mol/m2]
Flux(:,:) = ( Total_Spc / Total_Area ) * Mask(:,:,1)
Flux = Total_Spc / Total_Area

! Update species concentrations at surface [mol/mol]
State_Chm%Species(N)%Conc(:,:,1) = State_Chm%Species(N)%Conc(:,:,1) &
+ ( Flux(:,:) * AVO ) &
+ ( ( Flux * AVO ) &
/ ( State_Met%BXHEIGHT(:,:,1) &
* State_Met%AIRNUMDEN(:,:,1) )
* State_Met%AIRNUMDEN(:,:,1) ) ) * Mask(:,:,1)

ELSE
ErrMsg = TRIM( SpcInfo%Name ) // ': Src_Mode ' // &
Expand Down Expand Up @@ -532,17 +485,13 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &
REAL(fp) :: DT
REAL(fp) :: DecayConstant
REAL(fp) :: DecayRate
REAL(fp) :: LO3_kg

! Strings
CHARACTER(LEN=255) :: ErrMsg, ThisLoc

! Arrays
REAL(fp) :: Mask(State_Grid%NX,State_Grid%NY,State_Grid%NZ)

! Pointers to fields in the HEMCO data structure
REAL(f4), POINTER :: LOSS_O3(:,:,:) => NULL()

! Parameters
REAL(fp), PARAMETER :: ln2 = 0.693147181E+00_fp
REAL(fp), PARAMETER :: Day2Sec = 86400_fp
Expand All @@ -560,7 +509,6 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &
ThisLoc = &
' -> at Tracer_Sink_Phase (in module GeosCore/tracer_mod.F90)'
DT = GET_TS_CHEM()
Mask = 1.0_fp

!========================================================================
! Apply tracer sink
Expand All @@ -569,6 +517,9 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &
! Loop over species
DO N = 1, State_Chm%nAdvect

! Initialize mask for each species
Mask = 1.0_fp

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

Expand Down Expand Up @@ -608,34 +559,6 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &

ENDIF

!--------------------------------------------------------------------
! Get chemical loss (if needed)
!--------------------------------------------------------------------
IF ( TRIM(SpcInfo%Snk_Mode) == 'chemical_loss' ) THEN

! For stOX tracer get O3 loss [molec/cm3/s] from HEMCO
IF ( TRIM(SpcInfo%Name ) == 'stOX' ) THEN

CALL HCO_GetPtr( HcoState, 'O3_LOSS', LOSS_O3, RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Cannot get pointer to O3_LOSS!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF

ELSE

ErrMsg = 'Snk_Mode: chemical_loss not currently supported ' // &
'for ' // TRIM( SpcInfo%Name ) // &
'Please modify species_database.yml or add this ' // &
'capability in tracer_mod.F90.'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN

ENDIF

ENDIF

!---------------------------------------------------------------------
! Determine mask
!---------------------------------------------------------------------
Expand All @@ -649,7 +572,7 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &
DO I = 1, State_Grid%NX

! Set mask to zero outside of latitude zone
IF ( State_Grid%YMid(I,J) < SpcInfo%Snk_LatMin .and. &
IF ( State_Grid%YMid(I,J) < SpcInfo%Snk_LatMin .or. &
State_Grid%YMid(I,J) > SpcInfo%Snk_LatMax ) THEN
Mask(I,J,:) = 0.0_fp
ENDIF
Expand Down Expand Up @@ -739,47 +662,31 @@ SUBROUTINE Tracer_Sink_Phase( Input_Opt, State_Chm, State_Grid, &
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) * ( DecayRate * Mask(I,J,L) )
IF ( Mask(I,J,L) > 0 ) THEN
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) * DecayRate
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

ELSE IF ( TRIM(SpcInfo%Snk_Mode) == 'constant' ) THEN

WHERE( Mask > 0.0_fp )
State_Chm%Species(N)%Conc = SpcInfo%Snk_Value
ENDWHERE

ELSE IF ( TRIM(SpcInfo%Snk_Mode) == 'chemical_loss' ) THEN

! Special handling for stOX tracer
IF ( TRIM(SpcInfo%Name ) == 'stOX' ) THEN

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED )&
!$OMP PRIVATE( I, J, L, LO3_kg )&
!$OMP COLLAPSE( 3 )
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX

! Convert L(O3) from [molec/cm3/s] to [kg]
LO3_kg = LOSS_O3(I,J,L) & ! in molec/cm3/s
* (SpcInfo%MW_g * 1000.0_fp) / AVO & ! molec/cm3/s -> kg/m3/s
* State_Met%AIRVOL(I,J,L) & ! kg/m3/s -> kg/s
* DT ! kg/s -> kg

State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) - ( LO3_kg * MASK(I,J,L) )

ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

ENDIF
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED )&
!$OMP PRIVATE( I, J, L )&
!$OMP COLLAPSE( 3 )
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
IF ( Mask(I,J,L) > 0 ) THEN
State_Chm%Species(N)%Conc(I,J,L) = SpcInfo%Snk_Value
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

ENDIF

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -62,8 +62,6 @@ VerboseOnCores: root # Accepted values: root all
--> HTAP : false # 2008-2010
--> UNIFORM_CO : true
# ----- NON-EMISSIONS DATA ----------------------------------------------------
--> GLOBAL_O3 : true # 2010-2019
--> O3_LOSS : true # 2010-2019
--> OLSON_LANDMAP : true # 1985
--> YUAN_MODIS_LAI : true # 2000-2020
#------------------------------------------------------------------------------
Expand Down Expand Up @@ -196,20 +194,6 @@ VerboseOnCores: root # Accepted values: root all
* DELPDRY ./Restarts/GEOSChem.Restart.$YYYY$MM$DD_$HH$MNz.nc4 Met_DELPDRY $YYYY/$MM/$DD/$HH EY xyz 1 * - 1 1
)))GC_RESTART

#==============================================================================
# --- O3 concentrations ---
#==============================================================================
(((GLOBAL_O3
${RUNDIR_GLOBAL_O3}
)))GLOBAL_O3

#==============================================================================
# --- O3 loss rates ---
#==============================================================================
(((O3_LOSS
${RUNDIR_O3_LOSS}
)))O3_LOSS

#==============================================================================
# --- Olson land map masks ---
#==============================================================================
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ EXPID: ./OutputDir/GEOSChem
COLLECTIONS: 'Restart',
'RadioNuclide',
'SpeciesConc',
'Budget',
#'Budget',
'CloudConvFlux',
'DryDep',
'LevelEdgeDiags',
Expand Down