From 016ed0e59b15a845b20541f21fafd67f51a53317 Mon Sep 17 00:00:00 2001 From: Christopher Holmes Date: Thu, 19 Oct 2023 16:24:26 -0400 Subject: [PATCH 1/2] Capture Trop and PBL height changes in budgets Budget diagnostics now include mass tendency due to changes in * Tropopause height (TransportTrop budget) * PBL height (MixingPBL budget) The implementation could probably be improved as it currently assumes that Mixing occurs every dynamic time step. Signed-off-by: Christopher Holmes --- GeosCore/main.F90 | 117 +++++++++++++++++++++++++++++++++++++ Headers/state_diag_mod.F90 | 60 +++++++++++++++++++ 2 files changed, 177 insertions(+) diff --git a/GeosCore/main.F90 b/GeosCore/main.F90 index a0f72819d..83f5715f9 100644 --- a/GeosCore/main.F90 +++ b/GeosCore/main.F90 @@ -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, & @@ -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' ) @@ -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) @@ -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 @@ -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, & diff --git a/Headers/state_diag_mod.F90 b/Headers/state_diag_mod.F90 index d820981b6..e978a8def 100644 --- a/Headers/state_diag_mod.F90 +++ b/Headers/state_diag_mod.F90 @@ -80,10 +80,12 @@ MODULE State_Diag_Mod REAL(f8), POINTER :: BudgetEmisDryDepPBL (:,:,:) REAL(f8), POINTER :: BudgetTransportFull (:,:,:) REAL(f8), POINTER :: BudgetTransportTrop (:,:,:) + REAL(f8), POINTER :: BudgetTransportTropHeight(:,:,:) REAL(f8), POINTER :: BudgetTransportPBL (:,:,:) REAL(f8), POINTER :: BudgetMixingFull (:,:,:) REAL(f8), POINTER :: BudgetMixingTrop (:,:,:) REAL(f8), POINTER :: BudgetMixingPBL (:,:,:) + REAL(f8), POINTER :: BudgetMixingPBLHeight (:,:,:) REAL(f8), POINTER :: BudgetConvectionFull (:,:,:) REAL(f8), POINTER :: BudgetConvectionTrop (:,:,:) REAL(f8), POINTER :: BudgetConvectionPBL (:,:,:) @@ -774,10 +776,12 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%BudgetEmisDryDepPBL => NULL() State_Diag%BudgetTransportFull => NULL() State_Diag%BudgetTransportTrop => NULL() + State_Diag%BudgetTransportTropHeight => NULL() State_Diag%BudgetTransportPBL => NULL() State_Diag%BudgetMixingFull => NULL() State_Diag%BudgetMixingTrop => NULL() State_Diag%BudgetMixingPBL => NULL() + State_Diag%BudgetMixingPBLHeight => NULL() State_Diag%BudgetConvectionFull => NULL() State_Diag%BudgetConvectionTrop => NULL() State_Diag%BudgetConvectionPBL => NULL() @@ -1475,6 +1479,20 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%BudgetTransportTrop, & State_Chm, State_Diag, RC ) IF ( RC /= GC_SUCCESS ) RETURN + + arrayID = 'State_Diag%BudgetTransportTropHeight' + diagID = 'BudgetTransportTropHeight' + IF ( am_I_Root ) WRITE(6,20) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%BudgetTransportTropHeight( IM, JM, nAdvect ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%BudgetTransportTropHeight = 0.0_f8 + State_Diag%Archive_BudgetTransportTrop = .TRUE. + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%BudgetTransportTropHeight, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF ! PBL-only transport @@ -1535,6 +1553,20 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%BudgetMixingTrop, & State_Chm, State_Diag, RC ) IF ( RC /= GC_SUCCESS ) RETURN + + ! Second array needed for diagnosing changes in PBL height + arrayID = 'State_Diag%BudgetMixingPBLHeight' + diagID = 'BudgetMixingPBLHeight' + IF ( am_I_Root ) WRITE(6,20) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%BudgetMixingPBLHeight( IM, JM, nAdvect ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%BudgetMixingPBLHeight = 0.0_f8 + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%BudgetMixingPBLHeight, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF ! PBL-only mixing @@ -6430,6 +6462,13 @@ SUBROUTINE Cleanup_State_Diag( State_Diag, RC ) State_Diag%BudgetTransportTrop => NULL() ENDIF + IF ( ASSOCIATED( State_Diag%BudgetTransportTropHeight ) ) THEN + DEALLOCATE( State_Diag%BudgetTransportTropHeight, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%BudgetTransportTropHeight', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%BudgetTransportTropHeight => NULL() + ENDIF + IF ( ASSOCIATED( State_Diag%BudgetTransportPBL ) ) THEN DEALLOCATE( State_Diag%BudgetTransportPBL, STAT=RC ) CALL GC_CheckVar( 'State_Diag%BudgetTransportPBL', 2, RC ) @@ -6458,6 +6497,13 @@ SUBROUTINE Cleanup_State_Diag( State_Diag, RC ) State_Diag%BudgetMixingPBL => NULL() ENDIF + IF ( ASSOCIATED( State_Diag%BudgetMixingPBLHeight ) ) THEN + DEALLOCATE( State_Diag%BudgetMixingPBLHeight, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%BudgetMixingPBLHeight', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%BudgetMixingPBLHeight => NULL() + ENDIF + IF ( ASSOCIATED( State_Diag%BudgetConvectionFull ) ) THEN DEALLOCATE( State_Diag%BudgetConvectionFull, STAT=RC ) CALL GC_CheckVar( 'State_Diag%BudgetConvectionFull', 2, RC ) @@ -8063,6 +8109,13 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isRank ) Rank = 2 IF ( isTagged ) TagId = 'ADV' + ELSE IF ( TRIM( Name_AllCaps ) == 'BUDGETTRANSPORTTROPHEIGHT' ) THEN + IF ( isDesc ) Desc = 'Troposphere-only total mass rate of ' // & + 'change in tropopause height' + IF ( isUnits ) Units = 'kg s-1' + IF ( isRank ) Rank = 2 + IF ( isTagged ) TagId = 'ADV' + ELSE IF ( TRIM( Name_AllCaps ) == 'BUDGETTRANSPORTPBL' ) THEN IF ( isDesc ) Desc = 'PBL-only total mass rate of change ' // & ' in column for transport' @@ -8112,6 +8165,13 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isRank ) Rank = 2 IF ( isTagged ) TagId = 'ADV' + ELSE IF ( TRIM( Name_AllCaps ) == 'BUDGETMIXINGPBLHEIGHT' ) THEN + IF ( isDesc ) Desc = 'PBL-only total mass rate of change ' // & + ' in column for PBL height change' + IF ( isUnits ) Units = 'kg s-1' + IF ( isRank ) Rank = 2 + IF ( isTagged ) TagId = 'ADV' + ELSE IF ( TRIM( Name_AllCaps ) == 'BUDGETCONVECTIONFULL' ) THEN IF ( isDesc ) Desc = 'Total mass rate of change in column ' // & 'for convection' From c404489f7d10634608b8249b36d7931930e585b5 Mon Sep 17 00:00:00 2001 From: Christopher Holmes Date: Thu, 19 Oct 2023 15:33:52 -0400 Subject: [PATCH 2/2] Fix PBL level calculations Use local scale height and level thickness to determine PBL top level and pressure thickness. The code previously used a fixed global value for scale height, which is too large near poles and too small in tropics. This update lowers the PBL top level in the tropics and raises it near the poles. The change is typically 1 model level or less. Signed-off-by: Christopher Holmes --- GeosCore/pbl_mix_mod.F90 | 183 +++++++++++---------------------------- 1 file changed, 50 insertions(+), 133 deletions(-) diff --git a/GeosCore/pbl_mix_mod.F90 b/GeosCore/pbl_mix_mod.F90 index 3b780e2f2..5efa234fa 100644 --- a/GeosCore/pbl_mix_mod.F90 +++ b/GeosCore/pbl_mix_mod.F90 @@ -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 ! @@ -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) !================================================================= @@ -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 @@ -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 @@ -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