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

Extending Stratospheric Adjustment to GCClassic #2525

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Extending Stratospheric Adjustment to GCClassic #2525

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6b2df4c
Adding stratospheric adjustment capability to GCClassic.
cbarker211 Aug 7, 2024
bf19411
Floating point bug fix.
cbarker211 Aug 8, 2024
0608db4
Adding SEFDH.
cbarker211 Aug 9, 2024
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
106 changes: 106 additions & 0 deletions GeosCore/flexgrid_read_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ MODULE FlexGrid_Read_Mod
PUBLIC :: FlexGrid_Read_CN
PUBLIC :: FlexGrid_Read_A1
PUBLIC :: FlexGrid_Read_A3
PUBLIC :: FlexGrid_Read_Dyn
PUBLIC :: FlexGrid_Read_I3_1
PUBLIC :: FlexGrid_Read_I3_2
PUBLIC :: Copy_I3_Fields
Expand Down Expand Up @@ -522,6 +523,111 @@ END SUBROUTINE FlexGrid_Read_A1
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: FlexGrid_Read_Dyn
!
! !DESCRIPTION: Routine to read DynHeating variables and attributes from a NetCDF
! file containing 3-hr time-averaged (Dyn) data (based on flexgrid_read_a1)
!\\
!\\
! !INTERFACE:
!
SUBROUTINE FlexGrid_Read_Dyn( YYYYMMDD, HHMMSS, Input_Opt, State_Grid, &
State_Met )
!
! !USES:
!
USE Input_Opt_Mod, ONLY : OptInput
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
USE Get_Met_Mod
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: YYYYMMDD ! GMT date in YYYY/MM/DD format
INTEGER, INTENT(IN) :: HHMMSS ! GMT time in hh:mm:ss format
TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object
TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(MetState), INTENT(INOUT) :: State_Met ! Meteorology State object
!
! !REVISION HISTORY:
! 14 Mar 2024 - C. Barker - Initial version
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
! Scalars
INTEGER :: t_index ! Time index
CHARACTER(LEN=16) :: stamp ! Time and date stamp
CHARACTER(LEN=255) :: v_name ! netCDF variable name
CHARACTER(LEN=255) :: errMsg ! Error message
CHARACTER(LEN=255) :: caller ! Name of this routine

! Saved scalars
INTEGER, SAVE :: lastDate = -1 ! Stores last YYYYMMDD value
INTEGER, SAVE :: lastTime = -1 ! Stores last hhmmss value

! Arrays
REAL*4 :: Q(State_Grid%NX,State_Grid%NY,State_Grid%NZ) ! Temporary data array

!======================================================================
! Skip if we have already read data for this date & time
!======================================================================
IF ( YYYYMMDD == lastDate .and. HHMMSS == lastTime ) THEN
stamp = TimeStamp_String( YYYYMMDD, HHMMSS )
WRITE( 6, 20 ) stamp
20 FORMAT( ' - FLEXGRID DynHeating met fields for ', a, &
' have been read already' )
RETURN
ENDIF

!======================================================================
! Select the proper time slice
!======================================================================

! Name of this routine (for error printout)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cbarker211 can we remove these commented-out lines if they aren't needed?

!caller = "FlexGrid_Read_A3cld (flexgrid_read_mod.F90)"

! Stop w/ error if the t_index is not a multiple of 3.
!IF (MOD(HHMMSS, 030000) == 0) THEN
! WRITE( errMsg, 100 )
!100 FORMAT( 'Time_index value must be an integer!' )
! CALL ERROR_STOP( errMsg, caller )
!ENDIF

!======================================================================
! Get met fields from HEMCO
!======================================================================
! Read DynHeating
v_name = "DynHeating"
CALL Get_Met_3D( Input_Opt, State_Grid, Q, TRIM(v_name), t_index=1 )
State_Met%DynHeating = Q

! Echo info
stamp = TimeStamp_String( YYYYMMDD, HHMMSS )
WRITE( 6, 10 ) stamp
10 FORMAT( ' - Found all DynH met fields for ', a )

!======================================================================
! Cleanup and quit
!======================================================================

! Save date & time for next iteration
lastDate = YYYYMMDD
lastTime = HHMMSS

END SUBROUTINE FlexGrid_Read_Dyn
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: FlexGrid_Read_a3
!
! !DESCRIPTION: Convenience wrapper for the following routines which read
Expand Down
18 changes: 18 additions & 0 deletions GeosCore/hco_interface_gc_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4455,6 +4455,24 @@ SUBROUTINE Get_Met_Fields( Input_Opt, State_Chm, State_Grid, &
CALL FlexGrid_Read_A3( D(1), D(2), Input_Opt, State_Grid, State_Met )
ENDIF

!---------------------------------------------------------
! Read 3-hr time averaged DynHeating data (crb, 15/03/24)
!---------------------------------------------------------
IF (Input_Opt%RRTMG_FDH .and. Input_Opt%Read_Dyn_Heating) THEN
! We can use the boundary condition time functions
! as they also have no time adjustment
IF ( PHASE == 0 ) THEN
D = GET_FIRST_BC_TIME()
ELSE
D = GET_BC_TIME()
ENDIF
! We can use ITS_TIME_FOR_A3 as we are also loading in the data every 3 hrs.
IF ( PHASE == 0 .or. ITS_TIME_FOR_A3() .and. &
.not. ITS_TIME_FOR_EXIT() ) THEN
CALL FlexGrid_Read_Dyn( D(1), D(2), Input_Opt, State_Grid, State_Met )
ENDIF
ENDIF

!----------------------------------
! Read 3-hr instantanous data
!----------------------------------
Expand Down
7 changes: 0 additions & 7 deletions GeosCore/input_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2884,13 +2884,6 @@ SUBROUTINE Config_RRTMG( Config, Input_Opt, RC )
ENDIF
ENDIF

#ifndef MODEL_GCHPCTM
If (Input_Opt%RRTMG_FDH) Then
errMsg = 'Fixed dynamical heating in RRTMG is currently only available in GCHP'
CALL GC_Error( errMsg, RC, thisLoc )
End If
#endif

If (Input_Opt%RRTMG_SEFDH.and.(.not.Input_Opt%RRTMG_FDH)) Then
errMsg = 'Cannot have seasonally evolving FDH without enabling FDH!'
CALL GC_Error( errMsg, RC, thisLoc )
Expand Down
83 changes: 72 additions & 11 deletions Interfaces/GCClassic/main.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,6 +255,7 @@ PROGRAM GEOS_Chem
#ifdef RRTMG
! For stratospheric adjustment
REAL(f8), ALLOCATABLE :: DT_3D(:,:,:)
REAL(f8), ALLOCATABLE :: DT_3D_UPDATE(:,:,:)
REAL(f8), ALLOCATABLE :: HR_3D(:,:,:)
#endif

Expand Down Expand Up @@ -1658,6 +1659,42 @@ PROGRAM GEOS_Chem
WRITE( 6, '(a)' ) REPEAT( '#', 79 )
ENDIF

! Allocate temperature difference arrays (crb, 14/02/24)
If ( Input_Opt%RRTMG_FDH ) THEN
Allocate(DT_3D(State_Grid%NX,State_Grid%NY,State_Grid%NZ),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating DT_3D', ThisLoc )

! If using seasonally evolving FDH, need to grab the internal
! state temperature adjustment array
IF (Input_Opt%RRTMG_SEFDH) THEN
! DT_3D will be updated by the first call to RRTMG; also need to
! store the "current" value
Allocate(DT_3D_UPDATE(State_Grid%NX,State_Grid%NY,State_Grid%NZ),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating DT_3D_UPDATE', ThisLoc )
! Store the adjustment as previously projected to this time point
DT_3D(:,:,:) = State_Chm%TStrat_Adj(:,:,:)
! This will just hold the end-of-step value
DT_3D_UPDATE(:,:,:) = 0
ELSE
DT_3D(:,:,:) = 0
END IF

Allocate(HR_3D(State_Grid%NX,State_Grid%NY,State_Grid%NZ),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating HR_3D', ThisLoc )
HR_3D(:,:,:) = 0

! Read in dynamical heating rates if necessary
IF (Input_Opt%Read_Dyn_Heating) THEN
HR_3D(:,:,:) = State_Met%DynHeating(:,:,:)
ENDIF
ELSE
! Safer
Allocate(DT_3D(0,0,0),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating DT_3D', ThisLoc )
Allocate(HR_3D(0,0,0),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating HR_3D', ThisLoc )
ENDIF

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cbarker211 Could this be placed into a subroutine after the CONTAINS statement in main.F90 and just called from here? That would reduce clutter in an already-long main.F90 program.

State_Chm%RRTMG_iSeed = State_Chm%RRTMG_iSeed + 15

!------------------------------------------------------------------
Expand All @@ -1676,9 +1713,10 @@ PROGRAM GEOS_Chem
! Calculation for each of the potential output types
! See: wiki.geos-chem.org/Coupling_GEOS-Chem_with_RRTMG
!
! RRTMG outputs (scheduled in HISTORY.rc):
! 0-BA 1=O3 2=ME 3=SU 4=NI 5=AM
! 6=BC 7=OA 8=SS 9=DU 10=PM 11=ST
! 0=BASE and then...
! 1=O3 2=O3T 3=ME 4=H2O 5=CO2 6=CFC 7=N2O
! 8=SU 9=NI 10=AM 11=BC 12=OA 13=SS 14=DU
! 15=PM 16=ST
!
! State_Diag%RadOutInd(1) will ALWAYS correspond to BASE due
! to how it is populated from HISTORY.rc diaglist_mod.F90.
Expand All @@ -1695,13 +1733,9 @@ PROGRAM GEOS_Chem
! Generate mask for species in RT
CALL Set_SpecMask( State_Diag%RadOutInd(N), State_Chm )

! Dummy values (FDH not available in GC-Classic)
Allocate(DT_3D(0,0,0),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating DT_3D', ThisLoc )
Allocate(HR_3D(0,0,0),Stat=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error allocating HR_3D', ThisLoc )

! Compute radiative transfer for the given output
! If FDH is used, this step will be used to calculate DT and
! fill out DT_3D. The same will be true for HR_3D
CALL Do_RRTMG_Rad_Transfer( ThisDay = Day, &
ThisMonth = Month, &
iCld = State_Chm%RRTMG_iCld, &
Expand All @@ -1727,8 +1761,10 @@ PROGRAM GEOS_Chem

! Calculate for rest of outputs, if any
DO N = 2, State_Diag%nRadOut
! This time around, DT_3D is read in but not overwritten
WRITE( 6, 520 ) State_Diag%RadOutName(N), State_Diag%RadOutInd(N)
CALL Set_SpecMask( State_Diag%RadOutInd(N), State_Chm )
! This call will NOT update DT_3D, so we can just reuse the array
CALL Do_RRTMG_Rad_Transfer( ThisDay = Day, &
ThisMonth = Month, &
iCld = State_Chm%RRTMG_iCld, &
Expand Down Expand Up @@ -1763,8 +1799,33 @@ PROGRAM GEOS_Chem
CALL Debug_Msg( '### MAIN: a DO_RRTMG_RAD_TRANSFER' )
ENDIF

If (Allocated(DT_3D)) Deallocate(DT_3D)
If (Allocated(HR_3D)) Deallocate(HR_3D)
! Copy the adjustment back to DT_3D as calculated in the baseline calculation
IF (Input_Opt%RRTMG_SEFDH) THEN
DT_3D(:,:,:) = DT_3D_UPDATE(:,:,:)
END IF

! Store temperature change THEN heating rate from RRTMG in diagnostics (crb, 14/02/24)
If (Input_Opt%RRTMG_FDH) Then
IF (State_Diag%Archive_DynHeating) THEN
State_Diag%DynHeating(:,:,:) = HR_3D(:,:,:)
ENDIF
! NB: DT_3D is the temperature adjustment either after equilibration (pure FDH)
! or at the start of the NEXT radiation time step (SEFDH)
IF (State_Diag%Archive_DTRad ) THEN
State_Diag%DTRad(:,:,:) = DT_3D(:,:,:)
ENDIF
IF (Input_Opt%RRTMG_SEFDH) THEN
State_Chm%TStrat_Adj(:,:,:) = DT_3D(:,:,:)
END IF
ENDIF

RC = 0
IF (Allocated(DT_3D)) Deallocate(DT_3D, STAT=RC)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be more efficient to allocate all necessary arrays once in initialization and then to deallocate them at finalaization rather than on every timestep.

IF ( RC /= 0 ) Call Error_Stop( 'Error deallocating DT_3D', ThisLoc )
IF (Allocated(HR_3D)) Deallocate(HR_3D, STAT=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error deallocating HR_3D', ThisLoc )
IF (Allocated(DT_3D_UPDATE)) Deallocate(DT_3D_UPDATE, STAT=RC)
IF ( RC /= 0 ) Call Error_Stop( 'Error deallocating DT_3D_UPDATE', ThisLoc )

IF ( Input_Opt%useTimers ) THEN
CALL Timer_End( "RRTMG", RC )
Expand Down