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.

Sign up for GitHub

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

Jump to bottom

12.9.3+bugfix budget diagnostic #2002

Closed
wants to merge 2 commits into from
Closed

12.9.3+bugfix budget diagnostic #2002

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
016ed0e
Capture Trop and PBL height changes in budgets
cdholmes Oct 19, 2023
c404489
Fix PBL level calculations
cdholmes Oct 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
117 changes: 117 additions & 0 deletions GeosCore/main.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1170,6 +1170,25 @@ PROGRAM GEOS_Chem
! Do not compute any data in dry-run mode
IF ( notDryRun ) THEN

!------------------------------------------------------
! Tropopause height budget diagnostics - Part 1 of 2
!------------------------------------------------------
IF ( State_Diag%Archive_BudgetTransport ) THEN
! Get initial column masses
CALL Compute_Column_Mass( Input_Opt, &
State_Chm, State_Grid, State_Met, &
State_Chm%Map_Advect, &
.False., &
State_Diag%Archive_BudgetTransportTrop, &
.False., &
State_Diag%BudgetMass1, &
RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Tropopause height budget diagnostics error 1'
CALL GC_Error( ErrMsg, RC, ThisLoc )
ENDIF
ENDIF

! Interpolate I-3 fields to current dynamic timestep,
! based on their values at NSEC and NSEC+N_DYN
CALL Interp( NSECb, ELAPSED_TODAY, N_DYN, &
Expand Down Expand Up @@ -1232,6 +1251,40 @@ PROGRAM GEOS_Chem
ErrMsg = 'Error encountered in "Get_Cosine_SZA"!'
CALL Error_Stop( ErrMsg, ThisLoc )
ENDIF

!------------------------------------------------------
! Tropopause height budget diagnostics - Part 2 of 2
!------------------------------------------------------
IF ( State_Diag%Archive_BudgetTransport ) THEN

! Get final column masses and compute diagnostics
CALL Compute_Column_Mass( Input_Opt, &
State_Chm, State_Grid, State_Met, &
State_Chm%Map_Advect, &
.False., &
State_Diag%Archive_BudgetTransportTrop, &
.False., &
State_Diag%BudgetMass2, &
RC )

CALL Compute_Budget_Diagnostics( State_Grid, &
State_Chm%Map_Advect, &
real(Get_Ts_Dyn(), f8), &
.False., &
State_Diag%Archive_BudgetTransportTrop, &
.False., &
State_Diag%BudgetTransportFull, & ! This array shouldn't be changed, but can we avoid passing it?
State_Diag%BudgetTransportTropHeight,&
State_Diag%BudgetTransportPBL, & ! This array shouldn't be changed, but can we avoid passing it?
State_Diag%BudgetMass1, &
State_Diag%BudgetMass2, &
RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Tropopause height budget diagnostics error 2'
CALL GC_Error( ErrMsg, RC, ThisLoc )
ENDIF
ENDIF

ENDIF

IF ( prtDebug ) CALL Debug_Msg( '### MAIN: a INTERP, etc' )
Expand Down Expand Up @@ -1408,6 +1461,25 @@ PROGRAM GEOS_Chem
CALL Timer_Start( "Boundary layer mixing", RC )
ENDIF

!------------------------------------------------------
! PBL height budget diagnostics - Part 1 of 2
!------------------------------------------------------
IF ( State_Diag%Archive_BudgetMixing ) THEN
! Get initial column masses
CALL Compute_Column_Mass( Input_Opt, &
State_Chm, State_Grid, State_Met, &
State_Chm%Map_Advect, &
.False., &
.False., &
State_Diag%Archive_BudgetMixingPBL, &
State_Diag%BudgetMass1, &
RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Mixing budget diagnostics error 1 (Compute_PBL_Height)'
CALL GC_Error( ErrMsg, RC, ThisLoc )
ENDIF
ENDIF

! Move this call from the PBL mixing routines because the PBL
! height is used by drydep and some of the emissions routines.
! (ckeller, 3/5/15)
Expand All @@ -1419,6 +1491,39 @@ PROGRAM GEOS_Chem
CALL Error_Stop( ErrMsg, ThisLoc )
ENDIF

!------------------------------------------------------
! PBL height budget diagnostics - Part 2 of 2
!------------------------------------------------------
IF ( State_Diag%Archive_BudgetMixing ) THEN

! Get final column masses and compute diagnostics
CALL Compute_Column_Mass( Input_Opt, &
State_Chm, State_Grid, State_Met, &
State_Chm%Map_Advect, &
.False., &
.False., &
State_Diag%Archive_BudgetMixingPBL, &
State_Diag%BudgetMass2, &
RC )
CALL Compute_Budget_Diagnostics( State_Grid, &
State_Chm%Map_Advect, &
real(Get_Ts_Dyn(), f8), &
.False., &
.False., &
State_Diag%Archive_BudgetMixingPBL, &
State_Diag%BudgetMixingFull, & ! Can we avoid passing this array?
State_Diag%BudgetMixingTrop, & ! Can we avoid passing this array?
State_Diag%BudgetMixingPBLHeight, &
State_Diag%BudgetMass1, &
State_Diag%BudgetMass2, &
RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Mixing budget diagnostics error 2 (non-local mixing)'
CALL GC_Error( ErrMsg, RC, ThisLoc )
ENDIF
ENDIF


IF ( Input_Opt%useTimers ) THEN
CALL Timer_End( "Boundary layer mixing", RC )
ENDIF
Expand Down Expand Up @@ -1810,6 +1915,18 @@ PROGRAM GEOS_Chem
CALL Timer_Start( "=> History (netCDF diags)", RC )
ENDIF

! Sum the two transport budget diagnostics
IF ( State_Diag%Archive_BudgetTransportTrop ) THEN
State_Diag%BudgetTransportTrop = State_Diag%BudgetTransportTrop + &
State_Diag%BudgetTransportTropHeight
ENDIF

! Sum the two PBL mixing budget diagnostics
IF ( State_Diag%Archive_BudgetMixingPBL ) THEN
State_Diag%BudgetMixingPBL = State_Diag%BudgetMixingPBL + &
State_Diag%BudgetMixingPBLHeight
ENDIF

! Set State_Diag arrays that rely on state at end of timestep
CALL Set_Diagnostics_EndofTimestep( Input_Opt, State_Chm, &
State_Diag, State_Grid, &
Expand Down
183 changes: 50 additions & 133 deletions GeosCore/pbl_mix_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ SUBROUTINE COMPUTE_PBL_HEIGHT( State_Grid, State_Met, RC )
! !USES:
!
USE ErrCode_Mod
USE PhysConstants ! Scale height
USE PhysConstants ! g0, Rd
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
!
Expand All @@ -301,7 +301,7 @@ SUBROUTINE COMPUTE_PBL_HEIGHT( State_Grid, State_Met, RC )
!
LOGICAL :: Bad_Sum
INTEGER :: I, J, L, LTOP
REAL(fp) :: BLTOP, BLTHIK, DELP
REAL(fp) :: BLTOP, BLTHIK, DELP, Lower_Edge_Height
REAL(fp) :: P(0:State_Grid%NZ)

!=================================================================
Expand All @@ -312,158 +312,74 @@ SUBROUTINE COMPUTE_PBL_HEIGHT( State_Grid, State_Met, RC )
RC = GC_SUCCESS
Bad_Sum = .FALSE.
State_Met%InPbl = .FALSE.
State_Met%F_Under_PBLTop = 0d0
State_Met%F_of_PBL = 0d0

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, P, BLTOP, BLTHIK, LTOP, DELP )
!$OMP PRIVATE( I, J, L, P, BLTOP, BLTHIK, LTOP, DELP ) &
!$OMP Private( Lower_Edge_Height )
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX

!----------------------------------------------
! Define pressure edges:
! P(L-1) = P at bottom edge of box (I,J,L)
! P(L ) = P at top edge of box (I,J,L)
!----------------------------------------------
! PBL height above surface, m
State_Met%PBL_Top_m(I,J) = State_Met%PBLH(I,J)

! Pressure at level edges [hPa]
DO L = 0, State_Grid%NZ
P(L) = State_Met%PEDGE(I,J,L+1)
ENDDO

!----------------------------------------------
! Find PBL top and thickness [hPa]
!----------------------------------------------

! BLTOP = pressure at PBL top [hPa]
! Use barometric law since PBL is in [m]
BLTOP = P(0) * EXP( -State_Met%PBLH(I,J) / SCALE_HEIGHT )

! BLTHIK is PBL thickness [hPa]
BLTHIK = P(0) - BLTOP

!----------------------------------------------
! Find model level where BLTOP occurs
!----------------------------------------------

! Initialize
LTOP = 0

! Loop over levels
DO L = 1, State_Grid%NZ

! Exit when we get to the PBL top level
IF ( BLTOP > P(L) ) THEN
LTOP = L
EXIT
ELSE
State_Met%InPbl(I,J,L) = .TRUE.
ENDIF

ENDDO

!----------------------------------------------
! Define various related quantities
!----------------------------------------------

! IMIX(I,J) is the level where the PBL top occurs at (I,J)
! IMIX(I,J)-1 is the number of whole levels below the PBL top
IMIX(I,J) = LTOP
! Height of lower edge above surface, m
Lower_Edge_Height = 0d0

! Fraction of the IMIXth level underneath the PBL top
FPBL(I,J) = 1e+0_fp - ( BLTOP - P(LTOP) ) / &
( P(LTOP-1) - P(LTOP) )
! Find PBL top level (L) and pressure (hPa)
Do L=1, State_Grid%NZ

! PBL top [model layers]
State_Met%PBL_TOP_L(I,J) = FLOAT( IMIX(I,J) - 1 ) + FPBL(I,J)
If ( Lower_Edge_Height + State_Met%BXHEIGHT(I,J,L) > State_Met%PBLH(I,J) ) then

! PBL top [hPa]
State_Met%PBL_TOP_hPa(I,J) = BLTOP
! PBL top is in this level
LTOP = L

! Zero PBL top [m] -- compute below
State_Met%PBL_TOP_m(I,J) = 0e+0_fp
! Pressure at the PBL top altitude, hPa
! Use pressure lapse equation: p(PBLH) = p(z1) * exp( -(PBLH-z1) / Scale_Height )
! p(z1) = State_Met%PEDGE(I,J,L) = Pressure at the lower level edge
! PBLH - z1 = (State_Met%PBLH(I,J,L) - Lower_Edge_Height) = Height above the lower level edge
! Scale_Height = Rd * Tv / g0
State_Met%PBL_Top_hPa(I,J) = State_Met%PEdge(I,J,L) * &
EXP( -(State_Met%PBLH(I,J) - Lower_Edge_Height) * g0 / ( Rd * State_Met%Tv(I,J,L) ) )

! PBL thickness [hPa]
State_Met%PBL_THICK(I,J) = BLTHIK
! Fraction of PBL mass in layer L, will be normalized below
State_Met%F_of_PBL(I,J,L) = State_Met%PEdge(I,J,L) - State_Met%PBL_Top_hPa(I,J)

! Fraction of the grid cell mass under PBL top
State_Met%F_Under_PBLTop(I,J,L) = State_Met%F_of_PBL(I,J,L) / &
( State_Met%PEdge(I,J,L) - State_Met%PEdge(I,J,L+1) )

!==============================================================
! Loop up to edge of chemically-active grid
!==============================================================
DO L = 1, State_Grid%MaxChemLev
! Model level of PBL top (integer+fraction). The top is within level CEILING(PBL_Top_L)
State_Met%PBL_Top_L(I,J) = (LTOP-1) + State_Met%F_Under_PBLTop(I,J,L)

! Thickness of grid box (I,J,L) [hPa]
DELP = P(L-1) - P(L)
! PBL Thickness from surface to top, hPa
State_Met%PBL_Thick(I,J) = State_Met%PEdge(I,J,1) - State_Met%PBL_Top_hPa(I,J)

IF ( L < IMIX(I,J) ) THEN
!! Exit Do loop after we found PBL top level
Exit

!--------------------------------------------
! (I,J,L) lies completely below the PBL top
!--------------------------------------------
Else

! Fraction of grid box (I,J,L) w/in the PBL
State_Met%F_OF_PBL(I,J,L) = DELP / BLTHIK
! Grid cell fully within PBL
State_Met%inPBL(I,J,L) = .True.

! Fraction of grid box (I,J,L) underneath PBL top
State_Met%F_UNDER_PBLTOP(I,J,L) = 1e+0_fp
! Fraction of the grid cell mass under PBL top
State_Met%F_Under_PBLTop(I,J,L) = 1.0d0

! PBL height [m]
State_Met%PBL_TOP_m(I,J) = State_Met%PBL_TOP_m(I,J) + &
State_Met%BXHEIGHT(I,J,L)
! Fraction of PBL mass in layer L, will be normalized below
State_Met%F_of_PBL(I,J,L) = State_Met%PEdge(I,J,L) - State_Met%PEdge(I,J,L+1)

ELSE IF ( L == IMIX(I,J) ) THEN
! Update lower edge height, m
Lower_Edge_Height = Lower_Edge_Height + State_Met%BXHeight(I,J,L)

!--------------------------------------------
! (I,J,L) straddles the PBL top
!--------------------------------------------
EndIf

EndDo

! Fraction of grid box (I,J,L) w/in the PBL
State_Met%F_OF_PBL(I,J,L) = ( P(L-1) - BLTOP ) / BLTHIK

! Fraction of grid box (I,J,L) underneath PBL top
State_Met%F_UNDER_PBLTOP(I,J,L) = FPBL(I,J)

! PBL height [m]
State_Met%PBL_TOP_m(I,J) = State_Met%PBL_TOP_m(I,J) + &
( State_Met%BXHEIGHT(I,J,L) * &
FPBL(I,J) )

ELSE

!--------------------------------------------
! (I,J,L) lies completely above the PBL top
!--------------------------------------------

! Fraction of grid box (I,J,L) w/in the PBL
State_Met%F_OF_PBL(I,J,L) = 0e+0_fp

! Fraction of grid box (I,J,L) underneath PBL top
State_Met%F_UNDER_PBLTOP(I,J,L) = 0e+0_fp

ENDIF

!### Debug
!IF ( I==23 .and. J==34 .and. L < 6 ) THEN
! PRINT*, '###--------------------------------------'
! PRINT*, '### COMPUTE_PBL_HEIGHT'
! PRINT*, '### I, J, L : ', I, J, L
! PRINT*, '### P(L-1) : ', P(L-1)
! PRINT*, '### P(L) : ', P(L)
! PRINT*, '### F_OF_PBL : ', State_Met%F_OF_PBL(I,J,L)
! PRINT*, '### F_UNDER_TOP : ', &
! State_Met%F_UNDER_PBLTOP(I,J,L)
! PRINT*, '### IMIX : ', IMIX(I,J)
! PRINT*, '### FPBL : ', FPBL(I,J)
! PRINT*, '### PBL_TOP_hPa : ', State_Met%PBL_TOP_hPa(I,J)
! PRINT*, '### PBL_TOP_L : ', State_Met%PBL_TOP_L(I,J)
! PRINT*, '### DELP : ', DELP
! PRINT*, '### BLTHIK : ', BLTHIK
! PRINT*, '### BLTOP : ', BLTOP
! PRINT*, '### BXHEIGHT : ', State_Met%BXHEIGHT(I,J,L)
! PRINT*, '### PBL_TOP_m : ', State_Met%PBL_TOP_m(I,J)
! PRINT*, '### other way m : ', &
! P(0) * EXP( -State_Met%PBL_TOP_hPa(I,J) / SCALE_HEIGHT )
!ENDIF

ENDDO
! Fraction of PBL mass in layer L, now normalize to sum of 1
State_Met%F_of_PBL(I,J,:) = State_Met%F_of_PBL(I,J,:) / State_Met%PBL_Thick(I,J)

! Error check
IF ( ABS( SUM( State_Met%F_OF_PBL(I,J,:) ) - 1.e+0_fp) > 1.e-3_fp) THEN
Expand All @@ -472,6 +388,7 @@ SUBROUTINE COMPUTE_PBL_HEIGHT( State_Grid, State_Met, RC )
Bad_Sum = .TRUE.
!$OMP END CRITICAL
ENDIF

ENDDO
ENDDO
!$OMP END PARALLEL DO
Expand All @@ -484,7 +401,7 @@ SUBROUTINE COMPUTE_PBL_HEIGHT( State_Grid, State_Met, RC )
ENDIF

! Model level where PBL top occurs
State_Met%PBL_MAX_L = MAXVAL( IMIX )
State_Met%PBL_MAX_L = MAXVAL( CEILING( State_Met%PBL_Top_L ) )

END SUBROUTINE COMPUTE_PBL_HEIGHT
!EOC
Expand Down
Loading