diff --git a/drivers/ccpp/catchem_types.F90 b/drivers/ccpp/catchem_types.F90
new file mode 100644
index 00000000..eb010952
--- /dev/null
+++ b/drivers/ccpp/catchem_types.F90
@@ -0,0 +1,44 @@
+!> \file catchem_types.F90
+!> \brief Container module for CATCHEM model data structures
+!>
+!> \details
+!> Defines the core data structures used by CATCHEM chemistry model including
+!> arrays for meteorological states, chemical species concentrations, emissions,
+!> and diagnostic outputs. Provides the fundamental data container types needed
+!> for chemistry calculations.
+!>
+!> \author Barry Baker
+!>
+!> \date 11/2024
+!>
+!> \ingroup catchem_ccpp_group
+!!!>
+
+module catchem_types
+
+  use CATChem, only: MetStateType, ChemStateType, EmisStateType, DiagStateType
+
+  implicit none
+  private
+
+
+  !> Container type for CATCHEM model data
+  !!
+  !! Holds arrays for meteorological, chemical, emission and diagnostic states
+  !! along with the horizontal dimension parameter
+  !!
+  !! \ingroup catchem_ccpp_group
+  !!!>
+  type, public :: catchem_container_type
+    ! Array dimensions
+    integer :: im = 0  !> Horizontal dimension
+
+    ! State arrays
+    type(MetStateType),  allocatable :: met_state(:)    !> Meteorological state array
+    type(ChemStateType), allocatable :: chem_state(:)   !> Chemical state array
+    type(EmisStateType), allocatable :: emis_state(:)   !> Emission state array
+    type(DiagStateType), allocatable :: diag_state(:)   !> Diagnostic state array
+
+  end type catchem_container_type
+
+end module catchem_types
diff --git a/drivers/ccpp/catchem_wrapper_utils.F90 b/drivers/ccpp/catchem_wrapper_utils.F90
new file mode 100644
index 00000000..2bb431c2
--- /dev/null
+++ b/drivers/ccpp/catchem_wrapper_utils.F90
@@ -0,0 +1,394 @@
+!> \file catchem_wrapper_utils.F90
+!> \brief CATCHEM-CCPP interface utilities module
+!>
+!> \details
+!> Provides wrapper utilities for interfacing CATCHEM chemistry model with CCPP
+!> framework. Handles data transformation and management between host model and
+!> CATCHEM chemistry calculations.
+!>
+!> \author Barry Baker
+!>
+!> \date 11/2024
+!>
+!> \ingroup catchem_ccpp_group
+!!!>
+module catchem_wrapper_utils
+
+    use CATChem, only: ConfigType, MetStateType, ChemStateType, &
+                      EmisStateType, DiagStateType, cc_read_config, &
+                      cc_allocate_metstate, cc_allocate_chemstate, &
+                      cc_allocate_diagstate, cc_allocate_emisstate, &
+                      cc_checkdeallocate, cc_checkallocate, CC_SUCCESS
+    use catchem_types, only: catchem_container_type
+
+    implicit none
+    private
+
+    public :: read_catchem_config
+    public :: allocate_catchem_container
+    public :: transform_ccpp_to_catchem
+
+contains
+
+    !> Reads and initializes the CATChem configuration from a specified file
+    !!
+    !! This subroutine reads configuration settings and initializes the states
+    !! for meteorology, emissions, and chemistry using the CATChem interface.
+    !!
+    !! \param[inout] config         Configuration type containing CATChem settings
+    !! \param[inout] catchem_states Grid state container for CATChem components
+    !! \param[in]    config_file    Path to the CATChem configuration file
+    !! \param[out]   errflg         Error flag (0=success, non-zero=failure)
+    !! \param[out]   errmsg         Error message if errflg is non-zero
+    !!
+    !! \ingroup catchem_ccpp_group
+    !!!>
+    subroutine read_catchem_config(config, catchem_states, config_file, errflg, errmsg)
+        type(ConfigType), intent(inout) :: config
+        type(catchem_container_type), intent(inout) :: catchem_states
+        character(len=*), intent(in) :: config_file
+        integer, intent(out) :: errflg
+        character(len=*), intent(out) :: errmsg
+
+        call cc_read_config(config, &
+                          catchem_states%MetState(1), &
+                          catchem_states%EmisState(1), &
+                          catchem_states%ChemState(1), &
+                          errflg, &
+                          config_file)
+        if (errflg /= 0) then
+            errmsg = 'Error reading CATChem configuration file: '//trim(config_file)
+        end if
+    end subroutine read_catchem_config
+
+    !> Allocates memory for the entire catchem container
+    !!
+    !! \param[in]    config         CATChem configuration
+    !! \param[inout] catchem_states Container for all CATChem states
+    !! \param[in]    im            Number of horizontal points
+    !! \param[in]    kme           Number of vertical levels
+    !! \param[in]    nsoil         Number of soil levels
+    !! \param[out]   errflg        Error flag (0=success, non-zero=failure)
+    !! \param[out]   errmsg        Error message if errflg is non-zero
+    !!
+    !! \ingroup catchem_ccpp_group
+    !!!>
+    subroutine allocate_catchem_container(config, catchem_states, im, kme, nsoil, nLandTypes, errflg, errmsg)
+        type(ConfigType), intent(in) :: config
+        type(catchem_container_type), intent(inout) :: catchem_states
+        integer, intent(in) :: im
+        integer, intent(in) :: kme
+        integer, intent(in) :: nsoil
+        integer, intent(in) :: nLandTypes
+        integer, intent(out) :: errflg
+        character(len=*), intent(out) :: errmsg
+
+        integer :: i
+
+        ! Initialize
+        errflg = 0
+        errmsg = ''
+
+        ! Allocate State arrays in the container
+        if (.not. allocated(catchem_states%MetState)) then
+            allocate(catchem_states%MetState(im), stat=errflg)
+            if (check_allocation_error('catchem_states%MetState', errflg, errmsg)) return
+
+            ! Initialize each MetState using CC API
+            do i = 1, im
+                call CC_Allocate_MetState(catchem_states%MetState(i), errflg)
+                if (errflg /= CC_SUCCESS) then
+                    errmsg = 'Error in CC_Allocate_MetState'
+                    return
+                endif
+            end do
+        endif
+
+        if (.not. allocated(catchem_states%ChemState)) then
+            allocate(catchem_states%ChemState(im), stat=errflg)
+            if (check_allocation_error('catchem_states%ChemState', errflg, errmsg)) return
+
+            ! Initialize each ChemState using CC API
+            do i = 1, im
+                call CC_Allocate_ChemState(catchem_states%ChemState(i), config%Species_File, MetState(i), errflg)
+                if (errflg /= CC_SUCCESS) then
+                    errmsg = 'Error in CC_Allocate_ChemState'
+                    return
+                endif
+            end do
+        endif
+
+        if (.not. allocated(catchem_states%DiagState)) then
+            allocate(catchem_states%DiagState(im), stat=errflg)
+            if (check_allocation_error('catchem_states%DiagState', errflg, errmsg)) return
+
+            ! Initialize each DiagState using CC API
+            do i = 1, im
+                call CC_Allocate_DiagState(catchem_states%DiagState(i), config, MetState(i), errflg)
+                if (errflg /= CC_SUCCESS) then
+                    errmsg = 'Error in CC_Allocate_DiagState'
+                    return
+                endif
+            end do
+        endif
+
+        if (.not. allocated(catchem_states%EmisState)) then
+            allocate(catchem_states%EmisState(im), stat=errflg)
+            if (check_allocation_error('catchem_states%EmisState', errflg, errmsg)) return
+
+            ! Initialize each EmisState using CC API
+            do i = 1, im
+                call CC_Allocate_EmisState(catchem_states%EmisState(i), MetState(i), errflg)
+                if (errflg /= CC_SUCCESS) then
+                    errmsg = 'Error in CC_Allocate_EmisState'
+                    return
+                endif
+            end do
+        endif
+
+
+    end subroutine allocate_catchem_container
+
+    !> Transforms CCPP meteorological arrays into CATChem meteorological and chemistry states
+    !!
+    !! \param[in] im         Number of horizontal grid points
+    !! \param[in] kme        Number of vertical levels
+    !! \param[in] temp       3D temperature field (K)
+    !! \param[in] spechum    3D specific humidity field (kg/kg)
+    !! \param[in] pfull      3D full level pressure (Pa)
+    !! \param[in] phalf      3D half level pressure (Pa)
+    !! \param[in] u          3D zonal wind (m/s)
+    !! \param[in] v          3D meridional wind (m/s)
+    !! \param[in] delp       3D pressure thickness (Pa)
+    !! \param[in] zh         3D geopotential height (m)
+    !! \param[in] kh         3D vertical diffusivity (m2/s)
+    !! \param[in] prsl       3D layer mean pressure (Pa)
+    !! \param[in] prslk      3D Exner function
+    !! \param[in] u10m       10m u wind (m/s)
+    !! \param[in] v10m       10m v wind (m/s)
+    !! \param[in] tskin      Surface skin temperature (K)
+    !! \param[in] ps         Surface pressure (Pa)
+    !! \param[in] precip     Precipitation rate (kg/m2/s)
+    !! \param[in] slmsk      Land-sea mask (1=land,0=sea,2=ice)
+    !! \param[in] snowh      Snow depth (m)
+    !! \param[in] vegtype    Vegetation type
+    !! \param[in] soiltyp    Soil type
+    !! \param[in] hf         Sensible heat flux (W/m2)
+    !! \param[in] ust        Friction velocity (m/s)
+    !! \param[in] zpbl       PBL height (m)
+    !! \param[in] coszen     Cosine of solar zenith angle
+    !! \param[in] albedo     Surface albedo
+    !! \param[in] emis       Surface emissivity
+    !! \param[in] ustar      Friction velocity (m/s)
+    !! \param[in] shflx      Surface sensible heat flux (W/m2)
+    !! \param[in] lhflx      Surface latent heat flux (W/m2)
+    !! \param[in] snowc      Snow cover fraction (0-1)
+    !! \param[in] vegfrac    Green vegetation fraction (0-1)
+    !! \param[in] swdn       Surface downward SW radiation (W/m2)
+    !! \param[in] swup       Surface upward SW radiation (W/m2)
+    !! \param[in] lwdn       Surface downward LW radiation (W/m2)
+    !! \param[in] lwup       Surface upward LW radiation (W/m2)
+    !! \param[in] swdnc      Clear-sky surface downward SW radiation (W/m2)
+    !! \param[in] swupc      Clear-sky surface upward SW radiation (W/m2)
+    !! \param[in] lwdnc      Clear-sky surface downward LW radiation (W/m2)
+    !! \param[in] lwupc      Clear-sky surface upward LW radiation (W/m2)
+    !! \param[inout] MetState   CATChem meteorology state
+    !! \param[inout] ChemState  CATChem chemistry state
+    !! \param[out] errmsg     Error message
+    !! \param[out] errflg     Error flag
+    !!
+    !! \ingroup catchem_ccpp_group
+    !!!>
+    subroutine transform_ccpp_to_catchem(nVert, nVertInterface, nsoil, nlndcat, nsoilcat, lat, lon, &      ! Grid Information
+                                        ktau, dt, jdate, tile_num, area, xcosz, &  ! Grid Information
+                                        aero_rad_freq_opt, aero_feedback_opt, plmrise_freq_opt, & ! Model Options
+                                        lwi, dluse, soiltyp, vegtype_frac, &  ! Model Options
+                                        im, kme, nVert, nVertInterface, nsoil, nlndcat, nsoilcat, &  ! Model Options
+                                        temp, spechum, pfull, phalf, &  ! Meteorological Variables
+                                        u, v, delp, zh, kh, prsl, prslk, &  ! Meteorological Variables
+                                        u10m, v10m, tskin, ps, precip, &  ! Meteorological Variables
+                                        slmsk, snowh, vegtype, soiltyp, &  ! Surface Variables
+                                        hf, zpbl, coszen, albedo, emis, &  ! Surface Variables
+                                        ustar, shflx, lhflx, &  ! Near-Surface Meteorology
+                                        snowc, vegfrac, &  ! Surface Variables
+                                        swdn, swup, lwdn, lwup, &  ! Radiation Fluxes
+                                        swdnc, swupc, lwdnc, lwupc, &  ! Radiation Fluxes
+                                        MetState, ChemState, EmisState, DiagState, &  ! CATChem States
+                                        errmsg, errflg)  ! Error Handling
+                                        MetState, ChemState, EmisState, DiagState)
+
+      use CATChem, only: MetStateType, ChemStateType, DiagStateType, EmisStateType
+
+      implicit none
+      !! Transform CCPP meteorological arrays to CATChem states
+      integer,  intent(in)    :: im              !> number of horizontal points
+      integer,  intent(in)    :: kme             !> number of vertical levels
+      integer, intent(in)     :: nVert           !> number of vertical levels
+      integer, intent(in)     :: nVertInterface  !> number of vertical levels
+      integer, intent(in)     :: nsoil           !> number of soil levels
+      integer, intent(in)     :: nlndcat         !> number of land categories
+      integer, intent(in)     :: nsoilcat        !> number of soil categories
+
+      ! 3D/Layer Variables (dim(:,:))
+      real(kind=phys), intent(in) :: temp(:,:)       !> temperature (K)
+      real(kind=phys), intent(in) :: spechum(:,:)    !> specific humidity (kg/kg)
+      real(kind=phys), intent(in) :: pfull(:,:)      !> full level pressure (Pa)
+      real(kind=phys), intent(in) :: phalf(:,:)      !> half level pressure (Pa)
+      real(kind=phys), intent(in) :: u(:,:)          !> zonal wind (m/s)
+      real(kind=phys), intent(in) :: v(:,:)          !> meridional wind (m/s)
+      real(kind=phys), intent(in) :: delp(:,:)       !> pressure thickness (Pa)
+      real(kind=phys), intent(in) :: zh(:,:)         !> geopotential height (m)
+      real(kind=phys), intent(in) :: kh(:,:)         !> vertical diffusivity (m2/s)
+      real(kind=phys), intent(in) :: prsl(:,:)       !> layer mean pressure (Pa)
+      real(kind=phys), intent(in) :: prslk(:,:)      !> Exner function
+
+      ! Surface Variables (dim(:))
+      real(kind=phys), intent(in) :: ps(:)           !> surface pressure (Pa)
+      real(kind=phys), intent(in) :: tskin(:)        !> skin temperature (K)
+      real(kind=phys), intent(in) :: slmsk(:)        !> land-sea mask (1=land,0=sea,2=ice)
+      real(kind=phys), intent(in) :: snowh(:)        !> snow depth (m)
+      real(kind=phys), intent(in) :: vegtype(:)      !> vegetation type
+      real(kind=phys), intent(in) :: soiltyp(:)      !> soil type
+      real(kind=phys), intent(in) :: snowc(:)        !> snow cover fraction (0-1)
+      real(kind=phys), intent(in) :: vegfrac(:)      !> green vegetation fraction (0-1)
+
+      ! Near-Surface Meteorology (dim(:))
+      real(kind=phys), intent(in) :: u10m(:)         !> 10m u wind (m/s)
+      real(kind=phys), intent(in) :: v10m(:)         !> 10m v wind (m/s)
+      real(kind=phys), intent(in) :: ustar(:)        !> friction velocity (m/s)
+      real(kind=phys), intent(in) :: zpbl(:)         !> PBL height (m)
+
+      ! Surface Fluxes (dim(:))
+      real(kind=phys), intent(in) :: hf(:)           !> sensible heat flux (W/m2)
+      real(kind=phys), intent(in) :: shflx(:)        !> surface sensible heat flux (W/m2)
+      real(kind=phys), intent(in) :: lhflx(:)        !> surface latent heat flux (W/m2)
+      real(kind=phys), intent(in) :: precip(:)       !> precipitation rate (kg/m2/s)
+
+      ! Radiation Properties (dim(:))
+      real(kind=phys), intent(in) :: coszen(:)       !> cosine of solar zenith angle
+      real(kind=phys), intent(in) :: albedo(:)       !> surface albedo
+      real(kind=phys), intent(in) :: emis(:)         !> surface emissivity
+
+      ! Radiation Fluxes (dim(:))
+      real(kind=phys), intent(in) :: swdn(:)         !> downward shortwave radiation at surface (W/m2)
+      real(kind=phys), intent(in) :: swup(:)         !> upward shortwave radiation at surface (W/m2)
+      real(kind=phys), intent(in) :: lwdn(:)         !> downward longwave radiation at surface (W/m2)
+      real(kind=phys), intent(in) :: lwup(:)         !> upward longwave radiation at surface (W/m2)
+      real(kind=phys), intent(in) :: swdnc(:)        !> clear-sky downward shortwave radiation (W/m2)
+      real(kind=phys), intent(in) :: swupc(:)        !> clear-sky upward shortwave radiation (W/m2)
+      real(kind=phys), intent(in) :: lwdnc(:)        !> clear-sky downward longwave radiation (W/m2)2)
+      real(kind=phys), intent(in) :: lwupc(:)        !> clear-sky upward longwave radiation (W/m2)
+
+      ! Emissions
+      real(kind=phys), intent(in) :: emi_in(:)         !> emissions
+      ! CATChem States
+      type(MetStateType),  intent(inout) :: MetState(:)    !> CATChem meteorology state
+      type(ChemStateType), intent(inout) :: ChemState(:)   !> CATChem chemistry state
+      type(EmisStateType), intent(inout) :: EmisState(:)   !> CATChem emission state
+      type(DiagStateType), intent(inout) :: DiagState(:)   !> CATChem diagnostic state
+
+      ! Error handling
+      character(len=*), intent(out) :: errmsg    !> error message
+      integer,          intent(out) :: errflg    !> error flag
+
+      ! Local variables
+      integer :: i, k
+
+      ! Initialize error handling
+      errmsg = ''
+      errflg = 0
+
+      ! Check dimensions
+      if (size(MetState) /= im) then
+      errmsg = 'MetState dimension mismatch'
+      errflg = 1
+      return
+      endif
+
+      ! Transform data for each horizontal point
+      horiz: do i = 1, im
+        ! Verify vertical dimension
+        if (MetState(i)%nLEVS /= kme) then
+        errmsg = 'Vertical dimension mismatch'
+        errflg = 1
+        return
+        endif
+
+        ! 3D/Layer Variables
+        vert: do k = 1, kme
+            MetState(i)%temp(k)    = temp(i,k)
+            MetState(i)%spechum(k) = spechum(i,k)
+            MetState(i)%pfull(k)   = pfull(i,k)
+            MetState(i)%phalf(k)   = phalf(i,k)
+            MetState(i)%u(k)       = u(i,k)
+            MetState(i)%v(k)       = v(i,k)
+            MetState(i)%delp(k)    = delp(i,k)
+            MetState(i)%zh(k)      = zh(i,k)
+            MetState(i)%kh(k)      = kh(i,k)
+            MetState(i)%prsl(k)    = prsl(i,k)
+            MetState(i)%prslk(k)   = prslk(i,k)
+        end do vert
+
+        ! Surface Variables
+        MetState(i)%ps       = ps(i)
+        MetState(i)%tskin    = tskin(i)
+        MetState(i)%slmsk    = slmsk(i)
+        MetState(i)%snowh    = snowh(i)
+        MetState(i)%vegtype  = vegtype(i)
+        MetState(i)%soiltyp  = soiltyp(i)
+        MetState(i)%snowc    = snowc(i)
+        MetState(i)%vegfrac  = vegfrac(i)
+
+        ! Near-Surface Meteorology
+        MetState(i)%u10m     = u10m(i)
+        MetState(i)%v10m     = v10m(i)
+        MetState(i)%zpbl     = zpbl(i)
+        MetState(i)%ustar    = ustar(i)
+
+        ! Surface Fluxes
+        MetState(i)%hf       = hf(i)
+        MetState(i)%shflx    = shflx(i)
+        MetState(i)%lhflx    = lhflx(i)
+        MetState(i)%precip   = precip(i)
+        MetState(i)%ustar    = ustar(i)
+
+        ! Radiation Properties
+        MetState(i)%coszen   = coszen(i)
+        MetState(i)%albedo   = albedo(i)
+        MetState(i)%emis     = emis(i)
+
+        ! Radiation Fluxes
+        MetState(i)%swdn     = swdn(i)
+        MetState(i)%swup     = swup(i)
+        MetState(i)%lwdn     = lwdn(i)
+        MetState(i)%lwup     = lwup(i)
+        MetState(i)%swdnc    = swdnc(i)
+        MetState(i)%swupc    = swupc(i)
+        MetState(i)%lwdnc    = lwdnc(i)
+        MetState(i)%lwupc    = lwupc(i)
+
+      end do horiz
+
+    end subroutine transform_ccpp_to_catchem
+
+    !> Checks for allocation errors and sets error message
+    !!
+    !! \param[in] state_name Name of the state being allocated
+    !! \param[in] errflg     Error flag from allocation
+    !! \param[out] errmsg    Output error message if allocation failed
+    !! \return has_error     True if allocation error occurred
+    !! \ingroup catchem_ccpp_group
+    !!!>
+    function check_allocation_error(state_name, errflg, errmsg) result(has_error)
+        character(len=*), intent(in) :: state_name
+        integer, intent(in) :: errflg
+        character(len=*), intent(out) :: errmsg
+        logical :: has_error
+
+        has_error = (errflg /= 0)
+        if (has_error) then
+            errmsg = 'Error allocating '//trim(state_name)//' - catchem_wrapper_init'
+        end if
+    end function check_allocation_error
+
+end module catchem_wrapper_utils
diff --git a/drivers/ccpp/ccpp_catchem_interface.F90 b/drivers/ccpp/ccpp_catchem_interface.F90
new file mode 100644
index 00000000..2410b76d
--- /dev/null
+++ b/drivers/ccpp/ccpp_catchem_interface.F90
@@ -0,0 +1,774 @@
+!> \file catchem_interface_utils.F90
+!! \brief CATCHEM-CCPP interface utilities module
+!!
+!! \details
+!! This is the CCPP-Compliant wrapper for interfacing CATCHEM chemistry model with CCPP
+!! framework. Handles data transformation and management between host model and
+!! CATCHEM chemistry calculations. This module contains the init run and finalize functions and
+!! subroutines that facilitate the data exchange and coordinate transformations
+!! required for CATCHEM integration within CCPP-compliant models.
+!!
+!! \author Barry Baker, NOAA/OAR/ARL
+!!
+!! \date November 2024
+!! \defgroup catchem_ccpp_group CATChem CCPP Interface
+!! \ingroup catchem_ccpp_group
+!!
+!! \note This is part of the CATCHEM-CCPP interface layer that enables
+!!       chemistry calculations within the CCPP framework
+!!!>
+module catchem_interface
+
+  use catchem_types, only: catchem_container_type !> CATChem container type
+
+  use physcons, only: g => con_g, pi => con_pi
+  use machine, only: kind_phys
+  use catchem_config
+  use CATChem
+  use catchem_wrapper_utils
+
+implicit none
+
+private
+
+public :: catchem_interface_init, catchem_runphase1_interface_run
+
+type(ConfigType) :: Config                          !> CATChem configuration object
+type(catchem_container_type) :: CATChemStates       !> Container for all CATChem states
+
+!   integer :: im    !> Number of horizontal points
+!   integer :: kme   !> Number of vertical levels
+!   integer :: nsoil !> Number of soil layers
+
+contains
+
+
+
+   !> \section arg_table_catchem_interface_init Argument Table
+   !! \htmlinclude catchem_interface_init.html
+   !!
+   !! \brief Initialize the CATChem container
+   !! \param[in] catchem_configfile_in Name of the CATChem configuration file
+   !! \param[in] do_catchem_in Flag to enable CATChem calculations
+   !! \param[in] export_catchem_diags_in Flag to export CATChem diagnostics
+   !! \param[in] n_dbg_lines_in Number of debug lines
+   !! \param[in] im Number of horizontal points
+   !! \param[in] kme Number of vertical levels
+   !! \param[in] nsoil Number of soil layers
+   !! \param[out] errmsg Error message
+   !! \param[out] errflg Error flag
+   !!
+   !! \note This subroutine initializes the CATChem container and reads the configuration
+   !!       file to set up the chemistry model
+   !!
+   !! \ingroup catchem_ccpp_group
+   !!!>
+   subroutine catchem_interface_init(catchem_configfile_in, do_catchem, &
+                                 export_catchem_diags_in, n_dbg_lines_in, &
+                                 im, kme, nsoil, nLandType, errmsg, errflg)
+
+      implicit none
+
+      ! Input parameters
+      character(len=*), intent(in) :: catchem_configfile_in
+      logical,          intent(in) :: do_catchem
+      logical,          intent(in) :: export_catchem_diags_in
+      integer,          intent(in) :: n_dbg_lines_in
+      integer,          intent(in) :: im
+      integer,          intent(in) :: kme
+      integer,          intent(in) :: nsoil
+      integer,          intent(in) :: nLandType
+
+      ! Output parameters
+      character(len=*), intent(out) :: errmsg
+      integer,          intent(out) :: errflg
+
+      ! Local variables
+      integer :: i
+
+      errmsg = ''
+      errflg = 0
+
+      if (.not. do_catchem_in) return
+
+      ! Set global parameters
+      CATCHem_ConfigFile = catchem_config_opt
+      do_catchem = do_catchem_in
+      export_catchem_diags = export_catchem_diags_in
+      n_dbg_lines = n_dbg_lines_in
+
+      ! Initialize CATChem container
+      call CATChemStates%init(im)
+
+      ! Read configuration
+      call read_catchem_config(Config, CATChemStates, CATCHem_ConfigFile, errflg, errmsg)
+      if (errflg /= 0) return
+
+
+
+      ! Allocate states for each horizontal point
+      do i = 1, im
+         call allocate_states(Config, MetState(i), ChemState(i), DiagState(i), &
+                              EmisState(i), kme, nsoil, errflg, errmsg)
+         if (errflg /= 0) return
+      end do
+
+   end subroutine catchem_interface_init
+
+  !> \brief Brief description of the subroutine
+  !!
+  !! \section arg_table_catchem_gocart_interface_finalize Argument Table
+  !!
+  subroutine catchem_interface_finalize(im, kme, nsoil, errmsg, errflg)
+
+    use catchem_wrapper_utils, only: deallocate_states
+
+   implicit none
+
+   ! Input parameters
+   integer, intent(in) :: im
+   integer, intent(in) :: kme
+   integer, intent(in) :: nsoil
+
+   ! Output parameters
+   character(len=*), intent(out) :: errmsg
+   integer, intent(out) :: errflg
+
+   ! Local variables
+   integer :: i !> Horizontal index
+
+   do i = 1, im
+      call deallocate_states(MetState(i), ChemState(i), DiagState(i), EmisState(i), kme, nsoil, errflg, errmsg)
+      if (errflg /= 0) return
+   end do
+
+   deallocate(CATChemStates)
+
+  end subroutine catchem_wrapper_finalize
+
+  !>
+  !! This is the Configurable ATmospheric Chemistry (CATChem)
+  !! This is the CATChem interface Module
+  !! \section arg_table_catchem_runphase1_interface_run Argument Table
+  !! \htmlinclude catchem_runphase1_wrapper_run.html
+  !!
+  !>\section catchem_phase1_group CATChem Scheme General Algorithm
+  !> @{
+  subroutine catchem_interface_run(im, kte, kme, garea, nsoil, nlndcat, nsoilcat, &
+     rlat, rlon, tile_num, mpicomm, mpirank, mpiroot, & ! Grid information
+     do_catchem, & ! CATChem Flag on
+     ktau, dt, julian, jdate, idat, & ! Time information
+     xcosz, &
+     lwi, dluse, frlanduse, gvf, oro, frice, frocean, SoilMoist, SoilTemp & ! land water specific variables
+     ustar, u10m, v10m, tskin, t2m, dpt2m, hf2d, lf2d, znt, dswsfc, & ! land surface variables
+     pblh, kpbl, hpbl_thetav, soiltyp, tslb, snowdepth, frsnow, recmol, psim, psih, albedo, & ! land surface variables
+     psfc, prslp, pr3d, ph3d, phl3d, prl3d, tk3d, q3d, exch, us3d, vs3d, w, & ! State Variables
+     aer_ra_feedback_in, aer_ra_frq_in, & ! radiation
+     rainc_cpl, rain_cpl, cldf, chem_conv_tr_in, dqdt, & ! cloud variables
+     ntrac, ntchm, ntaero, chemarr_phys, chemarr, kemit_in, pert_scale_anthro, & ! Chemistry Variables
+     emis_amp_anthro, emi_in, dust_in, & ! Chemistry Variables
+     drydep, wetdpl, &  ! Chemistry Diagnostics
+     abem, pert_scale_plume, ca_emis_plume, biomass_burn_opt_in, plumerise_flag_in, & ! Fires
+     plumerisefire_freq_in, emis_amp_plume, ebu, fire_GBBEPx, fire_MODIS, & ! Fires
+     ca_sgs_emis, ca_sgs, & ! Cellular Automata
+     emis_amp_seas, do_sppt_emis, sppt_wts, ca_sgs_gbbepx_frp, emis_amp_dust, pert_scale_dust, & ! SPPT
+     errmsg, errflg)
+
+     implicit none
+
+     ! ARGUMENTS
+     !----------
+
+     ! CCPP Interface Variables
+     !-------------------------
+
+     ! Grid information
+     integer, intent(in) :: im       ! horizontal loop extent
+     integer, intent(in) :: kte      ! vertical layer dimension
+     integer, intent(in) :: kme      ! vertical interface dimension = number of vertical layers + 1
+     integer, intent(in) :: ktau     ! current forecast iteration
+     integer, intent(in) :: tile_num ! index of cube sphere tile
+     integer, intent(in) :: nsoil    !> vertical_dimension_of_soil
+     integer, intent(in) :: nlndcat !> number_of_vegetation_categories
+     integer, intent(in) :: nsoilcat !> number of soil categories
+     real(kind_phys), dimension(im), intent(in) :: garea !> grid area (m^2) of each grid cell
+     real(kind_phys), dimension(im), intent(in) :: rlat  !> latitude
+     real(kind_phys), dimension(im), intent(in) :: rlon, !> longitude
+     real(kind_phys), dimension(im), intent(in) :: xcosz !> cosine of solar zenith angle
+
+     ! Time information
+     integer, intent(in) :: idat(8)  ! initialization date and time (in iso order)
+     integer, intent(in) :: jdate    !> julian date
+     real(kind_phys), intent(in) :: dt      ! physics timestep (s)
+     real(kind_phys), intent(in) :: julian  ! forecast julian day (days)
+
+     ! MPI information
+     type(MPI_comm), intent(in) :: mpicomm
+     integer, intent(in) :: mpirank
+     integer, intent(in) :: mpiroot
+
+     ! Logicals
+     logical, intent(in) :: do_catchem
+
+     ! tracer information
+     integer, intent(in) :: ntrac    ! total number of tracers
+     integer, intent(in) :: ntchm      ! number of chemical tracers
+     integer, intent(in) :: ntaero      ! number of aerosol tracers
+     real(kind_phys), dimension(im, kte, ntrac), intent(inout) :: chemarr_phys
+     real(kind_phys), dimension(im, kte, ntrac), intent(inout) :: chemarr
+
+     ! Emissions
+     real(kind_phys), dimension(im, kte, 3), intent(in) :: emi_in ! 3 should be temporary... need to replace with either an input namelist option or have it be a pointer with variable dimension
+     real(kind_phys), dimension(im, kte, ntrac), intent(in) :: kemit_in !> emission inputs
+     real(kind_phys), dimension(im, kte), intent(in) :: dust_in         !> dust emission inputs
+     real(kind_phys), dimension(im), intent(in) :: pert_scale_anthro !> anthropogenic emission perturbation scale factor
+     real(kind_phys), dimension(im), intent(in) :: pert_scale_plume !> plume emission perturbation scale factor
+     real(kind_phys), dimension(im), intent(in) :: pert_scale_dust !> dust emission perturbation scale factor
+     real(kind_phys), dimension(im), intent(in) :: emis_amp_anthro !> anthropogenic emission amplitude
+     real(kind_phys), dimension(im), intent(in) :: emis_amp_plume !> plume emission amplitude
+     real(kind_phys), dimension(im), intent(in) :: emis_amp_seas !> seasonal emission amplitude
+     real(kind_phys), dimension(im), intent(in) :: emis_amp_dust !> dust emission amplitude
+     integer(kind_phys), intent(in) :: biomass_burn_opt_in !> biomass burning option
+     integer(kind_phys), intent(in) :: plumerise_flag_in !> plume rise flag
+     integer(kind_phys), intent(in) :: plumerisefire_freq_in !> plume rise fire frequency
+     real(kind_phys), dimension(im, kte), intent(inout) :: drydep !> dry deposition
+     real(kind_phys), dimension(im, kte), intent(inout) :: wetdpl !> wet deposition
+     real(kind_phys), dimension(im, kte), intent(inout) :: abem !> aerosol backscatter extinction mass
+     real(kind_phys), dimension(im, kte), intent(in) :: ca_emis_plume !> plume emission
+     logical, dimension(im, kte), intent(in) :: ca_sgs_emis !> SGS emission
+     logical, dimension(im, kte), intent(in) :: ca_sgs !> SGS
+     real(kind_phys), dimension(im, kte), intent(in) :: ca_sgs_gbbepx_frp !> SGS GBBEPx FRP
+     real(kind_phys), dimension(im, kte), intent(in) :: fire_GBBEPx !> GBBEPx fire emissions
+     real(kind_phys), dimension(im, kte), intent(in) :: fire_MODIS !> MODIS fire emissions
+     real(kind_phys), dimension(im, kte), intent(in) :: ebu !> EBU
+     real(kind_phys), dimension(im, kte), intent(in) :: sppt_wts !> SPPT weights
+     logical, intent(in) :: do_sppt_emis !> SPPT emission flag
+
+     ! land surface information
+     integer, dimension(im), intent(in)                :: lwi           !> sea land ice mask (sea = 0, land = 1, ice = 2)
+     integer, dimension(im), intent(in)                :: soiltyp       !> soil type
+     integer, dimension(im), intent(in)                :: dluse         !> vegetation type
+
+     real(kind_phys), dimension(im, nlndcat), intent(in) :: frlanduse     !> fraction of each land surface category
+     real(kind_phys), dimension(im), intent(in)        :: frice         !> fractin of ice cover
+     real(kind_phys), dimension(im), intent(in)        :: frocean       !> fraction of ocean cover
+     real(kind_phys), dimension(im), intent(in)        :: frsnow        !> fraction of snow cover over land
+     real(kind_phys), dimension(im), intent(in)        :: gvf           !> green vegetative fraction
+     real(kind_phys), dimension(im), intent(in)        :: oro           !> height above mean sea level (m)
+
+     real(kind_phys), dimension(im, nsoil), intent(in) :: soilmoist     !> volumetric fraction of soil moisture for lsm
+     real(kind_phys), dimension(im, nsoil), intent(in) :: soiltemp      !> soil temperature (K)
+     real(kind_phys), dimension(im,nsoil), intent(in)  :: tslb          !> soil temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: snowdepth     !> water equivalent snow depth (mm)
+     real(kind_phys), dimension(im), intent(in)        :: psfc          !> pressure at the surface (Pa)
+     real(kind_phys), dimension(im), intent(in)        :: prslp         !> sea level pressure (Pa)
+     real(kind_phys), dimension(im), intent(in)        :: pblh          !> PBL Thickness determined by the PBL scheme (m)
+     real(kind_phys), dimension(im), intent(in)        :: kpbl          !> PBL level
+     real(kind_phys), dimension(im), intent(in)        :: hpbl_thetav   !> PBL Height based on modified parcel method (m)
+     real(kind_phys), dimension(im), intent(in)        :: u10m          !> 10 m wind speed (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: v10m          !> 10 m wind speed (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: ustar         !> friction velocity (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: psim          !> Monin-Obukhov similarity parameter for momentum at 10m
+     real(kind_phys), dimension(im), intent(in)        :: psih          !> Monin-Obukhov similarity parameter for heat at 10m
+     real(kind_phys), dimension(im), intent(in)        :: tskin         !> skin temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: t2m           !> 2 m temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: dpt2m         !> 2 m dew point temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: hf2d          !> Sensible heat flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: lf2d          !> Latent heat flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: znt           !> surface roughness length in (cm)
+     real(kind_phys), dimension(im), intent(in)        :: dswsfc        !> downward short wave flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: recmol        !> one over obukhov length (m-1)
+     real(kind_phys), dimension(im), intent(in)        :: albedo        !> surface albedo
+
+     real(kind_phys), dimension(im, kte), intent(in) :: pr3d            !> air pressure at model layer interfaces (Pa)
+     real(kind_phys), dimension(im, kte), intent(in) :: prl3d           !> pressure at the model level (Pa)
+     real(kind_phys), dimension(im, kte), intent(in) :: ph3d            !> geopotential at the model level (m2 s-2)
+     real(kind_phys), dimension(im, kte), intent(in) :: phl3d           !> geopotential at the model layer interfaces (m2 s-2)
+     real(kind_phys), dimension(im, kte), intent(in) :: tk3d            !> temperature at the model level (K)
+     real(kind_phys), dimension(im, kte), intent(in) :: us3d            !> zonal wind at the model level (m/s)
+     real(kind_phys), dimension(im, kte), intent(in) :: vs3d            !> meridional wind at the model level (m/s)
+     real(kind_phys), dimension(im, kte), intent(in) :: q3d             !> specific humidity at the model level (kg/kg)
+     real(kind_phys), dimension(im, kte), intent(in) :: w               !> lagrangian_tendency_of_air_pressure
+     real(kind_phys), dimension(im, kte), intent(in) :: exch            !> atmospheric heat diffusivity
+
+     ! precipitation information
+     real(kind_phys), dimension(im), intent(in)        :: rain_cpl        !> total rain at this time step (m)
+     real(kind_phys), dimension(im), intent(in)        :: rainc_cpl       !> convective rain at this time step (m)
+     real(kind_phys), dimension(im), intent(in)        :: cldf            !> total cloud fraction
+     real(kind_phys), dimension(im, kte), intent(in)   :: dqdt            !> instantaneous_water_vapor_specific_humidity_tendency_due_to_convection
+     integer, intent(in)                               :: chem_conv_tr_in !> catchem convective transport option
+
+     ! Radiation
+     integer, intent(in) :: aer_ra_feedback_in  !> catchem aer radiation feedback option
+     integer, intent(in) :: aer_ra_frq_in       !> catchem_aer_ra_frq
+
+     ! Output
+     !-------
+     character(len=*), intent(out) :: errmsg
+     integer, intent(out) :: errflg
+
+     ! Local
+     !------
+     integer :: mpiid
+
+     ! Fill MetState and ChemState Arrays
+     call transform_ccpp_to_catchem(im, kte, kme, ntrac, ntc, ntr, &
+        gq0, qgrs, MetState, ChemState, DiagState, EmisState)
+
+     ! Run CATChem
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! -- volume to mass fraction conversion table (ppm -> ug/kg)
+     ppm2ugkg = 1._kind_phys
+     ppm2ugkg(p_sulf) = 1.e+03_kind_phys*mw_so4_aer/mwdry
+
+     ! -- set control flags
+     call_gocart = (mod(ktau, call_chemistry) == 0) .or. (ktau == 1)
+
+     ! -- compute accumulated large-scale and convective rainfall since last call
+     if (ktau > 1) then
+        dtstep = call_chemistry*dt
+     else
+        dtstep = dt
+     end if
+
+     !!!
+
+     !>- get ready for chemistry run
+     call catchem_gocart_prep( &
+        readrestart, chem_in_opt, ktau, dtstep, xcosz, &
+        garea, rlat, rlon, &
+        pr3d, ph3d, tk3d, prl3d, spechum, &
+        emi2_in, &
+        xlat, xlong, dxy, &
+        rri, t_phy, p_phy, rho_phy, dz8w, p8w, &
+        t8w, &
+        ntso2, ntsulf, ntDMS, ntmsa, ntpp25, &
+        ntbc1, ntbc2, ntoc1, ntoc2, ntpp10, &
+        ntrac, gq0, &
+        num_chem, num_moist, &
+        call_gocart, nvl_gocart, &
+        ttday, tcosz, gmt, julday, &
+        backg_oh, backg_h2o2, backg_no3, &
+        ppm2ugkg, &
+        moist, chem, &
+        ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte)
+     if (call_gocart) then
+        call gocart_chem_driver(ktau, dt, dtstep, gmt, julday, xcosz, &
+           t_phy, moist, chem, rho_phy, dz8w, p8w, backg_oh, oh_t, &
+           backg_h2o2, h2o2_t, backg_no3, no3_t, &
+           dxy, g, xlat, xlong, ttday, tcosz, &
+           chem_opt, num_chem, num_moist, &
+           ids, ide, jds, jde, kds, kde, &
+           ims, ime, jms, jme, kms, kme, &
+           its, ite, jts, jte, kts, kte)
+        call gocart_aerosols_driver(ktau, dtstep, t_phy, moist, &
+           chem, rho_phy, dz8w, p8w, dxy, g, &
+           chem_opt, num_chem, num_moist, &
+           ids, ide, jds, jde, kds, kde, &
+           ims, ime, jms, jme, kms, kme, &
+           its, ite, jts, jte, kts, kte)
+     end if
+
+     call sum_pm_gocart( &
+        rri, chem, pm2_5_dry, pm2_5_dry_ec, pm10, &
+        num_chem, chem_opt, &
+        ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte)
+
+     ! -- pm25 and pm10 for output , not for tracer options
+     do j = jts, jte
+        do k = kts, kte
+           do i = its, ite
+              pm25(i, j, k) = pm2_5_dry(i, k, j)
+              p10(i, j, k) = pm10(i, k, j)
+              !ebu_oc(i,j,k) = ebu      (i,k,j,p_ebu_oc)
+           end do
+        end do
+     end do
+
+     if (call_gocart) then
+        do j = jts, jte
+           do k = kts, kte
+              do i = its, ite
+                 oh_bg(i, j, k) = max(0., oh_t(i, k, j))
+                 h2o2_bg(i, j, k) = max(0., h2o2_t(i, k, j))
+                 no3_bg(i, j, k) = max(0., no3_t(i, k, j))
+              end do
+           end do
+        end do
+     end if
+
+     ! -- put chem stuff back into tracer array
+     do k = kts, kte
+        do i = its, ite
+           gq0(i, k, ntso2) = ppm2ugkg(p_so2)*max(epsilc, chem(i, k, 1, p_so2))
+           gq0(i, k, ntsulf) = ppm2ugkg(p_sulf)*max(epsilc, chem(i, k, 1, p_sulf))
+           gq0(i, k, ntdms) = ppm2ugkg(p_dms)*max(epsilc, chem(i, k, 1, p_dms))
+           gq0(i, k, ntmsa) = ppm2ugkg(p_msa)*max(epsilc, chem(i, k, 1, p_msa))
+           gq0(i, k, ntpp25) = ppm2ugkg(p_p25)*max(epsilc, chem(i, k, 1, p_p25))
+           gq0(i, k, ntbc1) = ppm2ugkg(p_bc1)*max(epsilc, chem(i, k, 1, p_bc1))
+           gq0(i, k, ntbc2) = ppm2ugkg(p_bc2)*max(epsilc, chem(i, k, 1, p_bc2))
+           gq0(i, k, ntoc1) = ppm2ugkg(p_oc1)*max(epsilc, chem(i, k, 1, p_oc1))
+           gq0(i, k, ntoc2) = ppm2ugkg(p_oc2)*max(epsilc, chem(i, k, 1, p_oc2))
+           gq0(i, k, ntpp10) = ppm2ugkg(p_p10)*max(epsilc, chem(i, k, 1, p_p10))
+        end do
+     end do
+
+     do k = kts, kte
+        do i = its, ite
+           qgrs(i, k, ntso2) = gq0(i, k, ntso2)
+           qgrs(i, k, ntsulf) = gq0(i, k, ntsulf)
+           qgrs(i, k, ntdms) = gq0(i, k, ntdms)
+           qgrs(i, k, ntmsa) = gq0(i, k, ntmsa)
+           qgrs(i, k, ntpp25) = gq0(i, k, ntpp25)
+           qgrs(i, k, ntbc1) = gq0(i, k, ntbc1)
+           qgrs(i, k, ntbc2) = gq0(i, k, ntbc2)
+           qgrs(i, k, ntoc1) = gq0(i, k, ntoc1)
+           qgrs(i, k, ntoc2) = gq0(i, k, ntoc2)
+           qgrs(i, k, ntpp10) = gq0(i, k, ntpp10)
+        end do
+     end do
+
+     !
+  end subroutine catchem_gocart_interface_run
+  !> @}
+
+  subroutine catchem_gocart_prep( &
+     readrestart, chem_in_opt, ktau, dtstep, xcosz, &
+     garea, rlat, rlon, &
+     pr3d, ph3d, tk3d, prl3d, spechum, &
+     emi2_in, &
+     xlat, xlong, dxy, &
+     rri, t_phy, p_phy, rho_phy, dz8w, p8w, &
+     t8w, &
+     ntso2, ntsulf, ntDMS, ntmsa, ntpp25, &
+     ntbc1, ntbc2, ntoc1, ntoc2, ntpp10, &
+     ntrac, gq0, &
+     num_chem, num_moist, &
+     call_gocart, nvl_gocart, &
+     ttday, tcosz, gmt, julday, &
+     backg_oh, backg_h2o2, backg_no3, &
+     ppm2ugkg, &
+     moist, chem, &
+     ids, ide, jds, jde, kds, kde, &
+     ims, ime, jms, jme, kms, kme, &
+     its, ite, jts, jte, kts, kte)
+
+     !Chem input configuration
+     logical, intent(in) :: readrestart
+     integer, intent(in) :: chem_in_opt, ktau, julday
+     real(kind=kind_phys), intent(in) :: dtstep, gmt
+
+     !FV3 input variables
+     integer, intent(in) :: ntrac
+     integer, intent(in) :: ntso2, ntpp25, ntbc1, ntoc1, ntpp10
+     integer, intent(in) :: ntsulf, ntbc2, ntoc2, ntDMS, ntmsa
+     real(kind=kind_phys), dimension(ims:ime), intent(in) :: garea, rlat, rlon, xcosz
+     real(kind=kind_phys), dimension(ims:ime, 64, 3), intent(in) :: emi2_in
+     real(kind=kind_phys), dimension(ims:ime, kms:kme), intent(in) :: pr3d, ph3d
+     real(kind=kind_phys), dimension(ims:ime, kts:kte), intent(in) :: &
+        tk3d, prl3d, spechum
+     real(kind=kind_phys), dimension(ims:ime, kts:kte, ntrac), intent(in) :: gq0
+
+     !GSD Chem variables
+     integer, intent(in) ::  num_chem, num_moist, &
+        nvl_gocart
+     logical, intent(in) ::  call_gocart
+     integer, intent(in) ::  ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte
+
+     real(kind_phys), dimension(num_chem), intent(in) :: ppm2ugkg
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: &
+        backg_oh, backg_h2o2, backg_no3
+
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: &
+        rri, t_phy, p_phy, rho_phy, dz8w, p8w, t8w
+     real(kind_phys), dimension(ims:ime, jms:jme), intent(out) :: &
+        xlat, xlong, dxy, &
+        ttday, tcosz
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_moist), intent(out) :: moist
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_chem), intent(out) :: chem
+
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme) :: z_at_w
+     real(kind_phys), dimension(nvl_gocart) :: p_gocart
+
+     ! -- local variables
+     real(kind_phys), dimension(ims:ime, jms:jme, nvl_gocart) :: oh_backgd, h2o2_backgd, no3_backgd
+     real(kind_phys) ::  pu, pl, aln, pwant
+     real(kind_phys) ::  xhour, xmin, gmtp, xlonn, xtime, real_time
+     real(kind_phys), DIMENSION(1, 1) :: sza, cosszax
+     integer i, ip, j, jp, k, kp, kk, kkp, nv, jmax, jmaxi, l, ll, n, ndystep, ixhour
+
+     p_gocart = (/1000., 992.5, 985., 977.5, 970., 955., 940., 925., 910., &
+        895., 880., 865., 850., 825., 800., 775., 750., 712.5, 675., 637.5, 600., &
+        562.5, 525., 487.5, 450., 412.5, 375., 337.5, 288.08, 244.88, 208.15, 176.93, &
+        150.39, 127.84, 108.66, 92.37, 78.51, 66.6, 56.39, 47.64, 40.18, 33.81, 28.37, &
+        23.73, 19.79, 16.46, 13.64, 11.28, 9.29, 7.62, 6.22, 5.05, 4.08, 3.28, 2.62, &
+        2.08, 1.65, 1.3, 1.02, 0.8, 0.62, 0.48, 0.37, 0.28/)
+
+     ! -- initialize output arrays
+     backg_oh = 0._kind_phys
+     backg_h2o2 = 0._kind_phys
+     backg_no3 = 0._kind_phys
+     rri = 0._kind_phys
+     t_phy = 0._kind_phys
+     p_phy = 0._kind_phys
+     rho_phy = 0._kind_phys
+     dz8w = 0._kind_phys
+     p8w = 0._kind_phys
+     t8w = 0._kind_phys
+     xlat = 0._kind_phys
+     xlong = 0._kind_phys
+     dxy = 0._kind_phys
+     ttday = 0._kind_phys
+     tcosz = 0._kind_phys
+     moist = 0._kind_phys
+     chem = 0._kind_phys
+     z_at_w = 0._kind_phys
+
+     do i = its, ite
+        dxy(i, 1) = garea(i)
+        xlat(i, 1) = rlat(i)*180./pi
+        xlong(i, 1) = rlon(i)*180./pi
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do i = its, ite
+           ip = i - its + 1
+           z_at_w(i, kts, j) = max(0., ph3d(ip, 1)/g)
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte
+           kp = k - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              dz8w(i, k, j) = abs(ph3d(ip, kp + 1) - ph3d(ip, kp))/g
+              z_at_w(i, k + 1, j) = z_at_w(i, k, j) + dz8w(i, k, j)
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte + 1
+           kp = k - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              p8w(i, k, j) = pr3d(ip, kp)
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte + 1
+           kk = min(k, kte)
+           kkp = kk - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              dz8w(i, k, j) = z_at_w(i, kk + 1, j) - z_at_w(i, kk, j)
+              t_phy(i, k, j) = tk3d(ip, kkp)
+              p_phy(i, k, j) = prl3d(ip, kkp)
+              rho_phy(i, k, j) = p_phy(i, k, j)/(287.04*t_phy(i, k, j)*(1.+.608*spechum(ip, kkp)))
+              rri(i, k, j) = 1./rho_phy(i, k, j)
+              moist(i, k, j, :) = 0.
+              moist(i, k, j, 1) = gq0(ip, kkp, p_atm_shum)
+              if (t_phy(i, k, j) > 265.) then
+                 moist(i, k, j, 2) = gq0(ip, kkp, p_atm_cldq)
+                 moist(i, k, j, 3) = 0.
+                 if (moist(i, k, j, 2) < 1.e-8) moist(i, k, j, 2) = 0.
+              else
+                 moist(i, k, j, 2) = 0.
+                 moist(i, k, j, 3) = gq0(ip, kkp, p_atm_cldq)
+                 if (moist(i, k, j, 3) < 1.e-8) moist(i, k, j, 3) = 0.
+              end if
+              !--
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        do k = 2, kte
+           do i = its, ite
+              t8w(i, k, j) = .5*(t_phy(i, k, j) + t_phy(i, k - 1, j))
+           end do
+        end do
+     end do
+
+     ! -- only used in phtolysis....
+     do j = jts, jte
+        do i = its, ite
+           t8w(i, 1, j) = t_phy(i, 1, j)
+           t8w(i, kte + 1, j) = t_phy(i, kte, j)
+        end do
+     end do
+
+     do k = kms, kte
+        do i = ims, ime
+           chem(i, k, jts, p_so2) = max(epsilc, gq0(i, k, ntso2)/ppm2ugkg(p_so2))
+           chem(i, k, jts, p_sulf) = max(epsilc, gq0(i, k, ntsulf)/ppm2ugkg(p_sulf))
+           chem(i, k, jts, p_dms) = max(epsilc, gq0(i, k, ntdms)/ppm2ugkg(p_dms))
+           chem(i, k, jts, p_msa) = max(epsilc, gq0(i, k, ntmsa)/ppm2ugkg(p_msa))
+           chem(i, k, jts, p_p25) = max(epsilc, gq0(i, k, ntpp25)/ppm2ugkg(p_p25))
+           chem(i, k, jts, p_bc1) = max(epsilc, gq0(i, k, ntbc1)/ppm2ugkg(p_bc1))
+           chem(i, k, jts, p_bc2) = max(epsilc, gq0(i, k, ntbc2)/ppm2ugkg(p_bc2))
+           chem(i, k, jts, p_oc1) = max(epsilc, gq0(i, k, ntoc1)/ppm2ugkg(p_oc1))
+           chem(i, k, jts, p_oc2) = max(epsilc, gq0(i, k, ntoc2)/ppm2ugkg(p_oc2))
+           chem(i, k, jts, p_p10) = max(epsilc, gq0(i, k, ntpp10)/ppm2ugkg(p_p10))
+        end do
+     end do
+
+     if (.NOT. readrestart) then
+        if (chem_in_opt == 0) then
+           if (ktau .le. 1) then
+              !           if(chem_opt > 0 ) then
+              do j = jts, jte
+                 jp = j - jts + 1
+                 do k = kts, kte
+                    do i = its, ite
+                       ip = i - its + 1
+                       if (chem_opt == CHEM_OPT_GOCART) then
+                          do n = 1, num_chem
+                             chem(i, k, j, n) = 1.e-30
+                          end do
+                       end if  ! chem_opt==300
+                       if ((chem_opt > CHEM_OPT_GOCART) .and. (chem_opt < CHEM_OPT_MAX)) then
+                          chem(i, k, j, p_so2) = 5.e-10
+                          chem(i, k, j, p_sulf) = 3.e-10
+                          chem(i, k, j, p_msa) = 1.e-10
+                          chem(i, k, j, p_dms) = 1.e-10
+                       end if !chem_opt >= 300 .and. chem_opt <  500
+
+                       !                if ((chem_opt == CHEM_OPT_GOCART_RACM) .or. (chem_opt == CHEM_OPT_RACM_SOA_VBS)) then  !added o3 background !lzhang
+                       !                  kk=min(k,kte)
+                       !                  kkp = kk - kts + 1
+                       !                  ! -- add initial constant into O3,CH4 and CO ect.
+                       !                  chem(i,k,j,p_o3)=epsilc
+                       !                  ! -- this section needs to be revisited before enabling the
+                       !                  ! corresponding chem_opt options
+                       !                  ! maxth=min(400.,th_pvsrf(i,j))
+                       !                  ! if (tr3d(ip,jp,kkp,p_atm_ptem) > maxth) then
+                       !                  !   chem(i,k,j,p_o3)=(airmw/48.)*tr3d(ip,jp,kkp,p_atm_o3mr)*1e6
+                       !                  !   !convert kg/kg to ppm
+                       !                  ! else
+                       !                  !   chem(i,k,j,p_o3)=0.03 !ppm
+                       !                  ! endif
+                       !                  chem(i,k,j,p_ch4)=1.85 !ppm
+                       !                  chem(i,k,j,p_co)=0.06 !ppm
+                       !                  chem(i,k,j,p_co2)=380.
+                       !                  chem(i,k,j,p_ete)=epsilc
+                       !                  chem(i,k,j,p_udd)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_hket)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_api)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_lim)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_dien)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_macr)=chem(i,k,j,p_ete)
+                       !                endif !( (chem_opt == 301.or.chem_opt==108))
+                    end do
+                 end do
+              end do
+           end if !(ktau<=1)
+
+        else !(chem_in_opt == 0 )
+
+           if ((ktau <= 1) .and. ((chem_opt == CHEM_OPT_GOCART_RACM) .or. (chem_opt == CHEM_OPT_RACM_SOA_VBS))) then  !added GFS o3 background above 380K!lzhang
+              do j = jts, jte
+                 jp = j - jts + 1
+                 do k = kts, kte + 1
+                    kk = min(k, kte)
+                    kkp = kk - kts + 1
+                    do i = its, ite
+                       ip = i - its + 1
+                       ! -- this section needs to be revisited before enabling the
+                       ! corresponding chem_opt options
+                       ! maxth=min(400.,th_pvsrf(i,j))
+                       ! if (tr3d(ip,jp,kkp,p_atm_ptem) >= maxth) then
+                       !   chem(i,k,j,p_o3)=(airmw/48.)*tr3d(ip,jp,kkp,p_atm_o3mr)*1e6 !convert kg/kg to ppm
+                       ! endif !380K
+                    end do
+                 end do
+              end do
+           end if ! chem_opt == 301.or.chem_opt==108
+
+        end if !(chem_in_opt == 1 )
+     end if ! readrestart
+
+     !-- assgin read in 3D background chemical species
+     do i = its, ite
+        do k = 1, nvl_gocart
+           h2o2_backgd(i, 1, k) = emi2_in(i, k, 1)
+           no3_backgd(i, 1, k) = emi2_in(i, k, 2)
+           oh_backgd(i, 1, k) = emi2_in(i, k, 3)
+        end do
+     end do
+
+     !
+     ! -- gocart background fields only if gocart is called
+     if (call_gocart .and. (chem_opt == CHEM_OPT_GOCART)) then
+        do j = jts, jte
+           do i = its, ite
+              do k = kts, kte
+                 do ll = 2, nvl_gocart
+                    l = ll
+                    if (p_gocart(l) < .01*p_phy(i, k, j)) exit
+                 end do
+                 pu = alog(p_gocart(l))
+                 pl = alog(p_gocart(l - 1))
+                 pwant = alog(.01*p_phy(i, k, j))
+                 if (pwant > pl) then
+                    backg_oh(i, k, j) = oh_backgd(i, j, l)
+                    backg_h2o2(i, k, j) = h2o2_backgd(i, j, l)
+                    backg_no3(i, k, j) = no3_backgd(i, j, l)
+                 else
+                    aln = (oh_backgd(i, j, l)*(pwant - pl) + &
+                       oh_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_oh(i, k, j) = aln
+                    aln = (h2o2_backgd(i, j, l)*(pwant - pl) + &
+                       h2o2_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_h2o2(i, k, j) = aln
+                    aln = (no3_backgd(i, j, l)*(pwant - pl) + &
+                       no3_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_no3(i, k, j) = aln
+                 end if
+              end do
+           end do
+        end do
+     end if   ! end gocart stuff
+     !endif !restart
+
+     if ((chem_opt == CHEM_OPT_RACM_SOA_VBS) .or. (chem_opt >= CHEM_OPT_GOCART .and. chem_opt < CHEM_OPT_MAX)) then
+        !ndystep=86400/ifix(dtstepc)
+        ndystep = 86400/ifix(dtstep)
+        do j = jts, jte
+           do i = its, ite
+              tcosz(i, j) = 0.
+              ttday(i, j) = 0.
+              xlonn = xlong(i, j)
+              do n = 1, ndystep
+                 xtime = n*dtstep/60.
+                 ixhour = ifix(gmt + .01) + ifix(xtime/60.)
+                 xhour = float(ixhour)
+                 xmin = 60.*gmt + (xtime - xhour*60.)
+                 gmtp = mod(xhour, 24.)
+                 gmtp = gmtp + xmin/60.
+                 CALL szangle(1, 1, julday, gmtp, sza, cosszax, xlonn, rlat(i))
+                 TCOSZ(i, j) = TCOSZ(I, J) + cosszax(1, 1)
+                 if (cosszax(1, 1) > 0.) ttday(i, j) = ttday(i, j) + dtstep
+              end do
+           end do
+        end do
+     end if !chem_opt >= 300 .and. chem_opt <  500
+
+  end subroutine catchem_gocart_prep
+
+  !> @}
+end module catchem_gocart_wrapper
diff --git a/drivers/ccpp/ccpp_catchem_interface.meta b/drivers/ccpp/ccpp_catchem_interface.meta
new file mode 100644
index 00000000..62a98298
--- /dev/null
+++ b/drivers/ccpp/ccpp_catchem_interface.meta
@@ -0,0 +1,758 @@
+[ccpp-table-properties]
+  name = catchem_interface
+  type = scheme
+  dependencies = ../../parameters/catchem_config.F90,../../parameters/catchem_constants.F90,../../src/dep_vertmx_mod.F90,../../src/drydep_gocart_mod.F90,../../src/drydep_wesely_mod.F90
+
+########################################################################
+[ccpp-arg-table]
+  name = catchem_interface_init
+  type = scheme
+[im]
+  standard_name = horizontal_loop_extent
+  long_name = horizontal loop extent
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[kte]
+  standard_name = vertical_layer_dimension
+  long_name = vertical layer dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[kme]
+  standard_name = vertical_interface_dimension
+  long_name = number of vertical levels plus one
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[nsoil]
+  standard_name = vertical_dimension_of_soil
+  long_name = soil vertical layer dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[nlndcat]
+  standard_name = number_of_vegetation_categories
+  long_name = number of vegetation categories
+  units = count
+  dimensions = ()
+  type = integer
+[nsoilcat]
+  standard_name = number_of_soil_categories
+  long_name = number of soil categories
+  units = count
+  dimensions = ()
+  type = integer
+[do_catchem]
+  standard_name = do_catchem_coupling
+  long_name = flag controlling catchem collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+########################################################################
+[ccpp-arg-table]
+  name = catchem_interface_finalize
+  type = scheme
+[im]
+  standard_name = horizontal_loop_extent
+  long_name = horizontal loop extent
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+
+########################################################################
+[ccpp-arg-table]
+  name = catchem_interface_run
+  type = scheme
+#Grid State
+[im]
+  standard_name = horizontal_loop_extent
+  long_name = horizontal loop extent
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[kte]
+  standard_name = vertical_layer_dimension
+  long_name = vertical layer dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[kme]
+  standard_name = vertical_interface_dimension
+  long_name = number of vertical levels plus one
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[ktau]
+  standard_name = index_of_timestep
+  long_name = current forecast iteration
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[garea]
+  standard_name = cell_area
+  long_name = grid cell area
+  units = m2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[nsoil]
+  standard_name = vertical_dimension_of_soil
+  long_name = soil vertical layer dimension
+  units = count
+  dimensions = ()
+  type = integer
+[nlndtype]
+  standard_name = number_of_vegetation_categories
+  long_name = number of vegetation categories
+  units = count
+  dimensions = ()
+  type = integer
+[nsoilcat]
+  standard_name = number_of_soil_categories
+  long_name = number of soil categories
+  units = count
+  dimensions = ()
+  type = integer
+[rlat]
+  standard_name = latitude
+  long_name = latitude
+  units = radian
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[rlon]
+  standard_name = longitude
+  long_name = longitude
+  units = radian
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[errmsg]
+  standard_name = ccpp_error_message
+  long_name = error message for error handling in CCPP
+  units = none
+  dimensions = ()
+  type = character
+  kind = len=*
+  intent = out
+[errflg]
+  standard_name = ccpp_error_code
+  long_name = error code for error handling in CCPP
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = out
+  # LOGICALS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[do_catchem]
+  standard_name = do_catchem_coupling
+  long_name = flag controlling catchem collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+# TIME VARIABLES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[dt]
+  standard_name = timestep_for_physics
+  long_name = physics time step
+  units = s
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[julian]
+  standard_name = forecast_julian_day
+  long_name = julian day
+  units = days
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[jdate]
+  standard_name = date_and_time_of_forecast_in_united_states_order
+  long_name = current forecast date and time
+  units = none
+  dimensions = (8)
+  type = integer
+  intent = in
+[idat]
+  standard_name = date_and_time_at_model_initialization_in_iso_order
+  long_name = initialization date and time
+  units = none
+  dimensions = (8)
+  type = integer
+  intent = in
+# Land Parameters #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[lwi]
+  standard_name = sea_land_ice_mask
+  long_name = landmask: sea/land/ice=0/1/2
+  units = flag
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = in
+[dluse]
+  standard_name = vegetation_type_classification
+  long_name = vegetation type at each grid cell
+  units = index
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = in
+[frlanduse]
+  standard_name = vegetation_type_classification
+  long_name = vegetation type at each grid cell
+  units = index
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = in
+[gvf]
+  standard_name = bounded_vegetation_area_fraction
+  long_name = areal fractional cover of green vegetation bounded on the bottom
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[fice]
+  standard_name = sea_ice_area_fraction_of_sea_area_fraction
+  long_name = ice fraction over open water
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[oceanfrac]
+  standard_name = sea_area_fraction
+  long_name = fraction of horizontal grid area occupied by ocean
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[lakefrac]
+  standard_name = flag_nonzero_lake_surface_fraction
+  long_name = flag indicating presence of some lake surface area fraction
+  units = flag
+  dimensions = (horizontal_loop_extent)
+  type = logical
+# Surface Met <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[ustar]
+  standard_name = surface_friction_velocity
+  long_name = boundary layer parameter
+  units = m s-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[u10m]
+  standard_name = x_wind_at_10m
+  long_name = 10 meter u wind speed
+  units = m s-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[v10m]
+  standard_name = y_wind_at_10m
+  long_name = 10 meter v wind speed
+  units = m s-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[tskin]
+  standard_name = surface_skin_temperature
+  long_name = surface skin temperature
+  units = K
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[hf2d]
+  standard_name = instantaneous_surface_upward_sensible_heat_flux
+  long_name = surface upward sensible heat flux valid for current call
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[pblh]
+  standard_name = atmosphere_boundary_layer_thickness
+  long_name = PBL thickness
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[dswsfc]
+  standard_name = surface_downwelling_shortwave_flux
+  long_name = surface downwelling shortwave flux at current time
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[znt]
+  standard_name = surface_roughness_length
+  long_name = surface roughness length in cm
+  units = cm
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[snowdepth]
+  standard_name = lwe_thickness_of_snow_amount_for_chemcoupling
+  long_name = total snow precipitation chem
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[frsnow]
+  standard_name = surface_snow_area_fraction_over_land
+  long_name = Surface snow area fraction over land
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[psih]
+  standard_name = Monin_Obukhov_similarity_function_for_heat_at_10m
+  long_name = Monin-Obukhov similarity parameter for heat at 10m
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[psim]
+  standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m
+  long_name = Monin-Obukhov similarity parameter for momentum at 10m
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[soilmoist]
+  standard_name = soil_temperature
+  long_name = volumetric fraction of soil moisture for lsm
+  units =
+  dimensions = (horizontal_loop_extent,vertical_dimension_of_soil)
+  type = real
+  kind = kind_phys
+  intent = inout
+[soiltemp]
+  standard_name = volume_fraction_of_condensed_water_in_soil
+  long_name = volumetric fraction of soil moisture for lsm
+  units = 'K'
+  dimensions = (horizontal_loop_extent,vertical_dimension_of_soil)
+  type = real
+  kind = kind_phys
+  intent = inout
+# State Variables <<<<<<<<<<<<<<<<<<<<<<<<<<<
+[prsfc]
+  standard_name = air_pressure_at_surface
+  long_name = Air pressure at local surface
+  units = Pa
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[prslp]
+  standard_name = us_standard_air_pressure_at_sea_level
+  long_name = standard atmospheric pressure at sea level
+  units = Pa
+  type = real
+  kind = kind_phys
+  intent = in
+[pr3d]
+  standard_name = air_pressure_at_interface
+  long_name = air pressure at model layer interfaces
+  units = Pa
+  dimensions = (horizontal_loop_extent,vertical_interface_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[ph3d]
+  standard_name = geopotential_at_interface
+  long_name = geopotential at model layer interfaces
+  units = m2 s-2
+  dimensions = (horizontal_loop_extent,vertical_interface_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[phl3d]
+  standard_name = geopotential
+  long_name = geopotential at model layer centers
+  units = m2 s-2
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[prl3d]
+  standard_name = air_pressure
+  long_name = mean layer pressure
+  units = Pa
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[tk3d]
+  standard_name = air_temperature_of_new_state
+  long_name = updated temperature
+  units = K
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[q3d]
+  standard_name = specific_humidity_of_new_state
+  long_name = water vapor specific humidity updated by physics
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[exch]
+  standard_name = atmosphere_heat_diffusivity
+  long_name = diffusivity for heat
+  units = m2 s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[us3d]
+  standard_name = x_wind_of_new_state
+  long_name = zonal wind updated by physics
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[vs3d]
+  standard_name = y_wind_of_new_state
+  long_name = meridional wind updated by physics
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[w]
+  standard_name = lagrangian_tendency_of_air_pressure
+  long_name = layer mean vertical velocity
+  units = Pa s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+# Radiation <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[aer_ra_feedback_in]
+  standard_name = catchem_aer_ra_feedback
+  long_name = catchem aer radiation feedback option
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[aer_ra_frq_in]
+  standard_name = catchem_aer_ra_frq
+  long_name = catchem aer radiation frequency
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+# Precipitation <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[rainc_cpl]
+  standard_name = cumulative_lwe_thickness_of_convective_precipitation_amount_for_coupling
+  long_name = total convective precipitation
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[rain_cplchm]
+  standard_name = lwe_thickness_of_precipitation_amount_for_chemcoupling
+  long_name = total rain precipitation chem
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[cldf]
+  standard_name = cloud_area_fraction
+  long_name = fraction of grid box area in which updrafts occur
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[chem_conv_tr_in]
+  standard_name = catchem_conv_tr
+  long_name = catchem convective transport option
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[dqdt]
+  standard_name = instantaneous_water_vapor_specific_humidity_tendency_due_to_convection
+  long_name = instantaneous moisture tendency due to convection
+  units = kg kg-1 s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+# Tracer info <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[ntrac]
+  standard_name = number_of_tracers
+  long_name = number of tracers
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[ntchm]
+  standard_name = number_of_chemical_tracers
+  long_name = number of chemical tracers
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[ntaero]
+  standard_name = number_of_aerosol_tracers
+  long_name = number of aerosol tracers
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[ntchmdiag]
+  standard_name = number_of_chemical_tracers_for_diagnostics
+  long_name = number of chemical tracers for diagnostic output
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[chemarr_phys]
+  standard_name = tracer_concentration_of_new_state
+  long_name = tracer concentration updated by physics
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers)
+  type = real
+  kind = kind_phys
+  intent = inout
+[chemarr]
+  standard_name = tracer_concentration
+  long_name = model layer mean tracer concentration
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers)
+  type = real
+  kind = kind_phys
+  intent = inout
+# Fires <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[abem]
+  standard_name = instantaneous_anthopogenic_and_biomass_burning_emissions
+  long_name = instantaneous anthopogenic and biomass burning emissions for black carbon, organic carbon, and sulfur dioxide
+  units = kg m-2 s-1
+  dimensions = (horizontal_loop_extent,12)
+  type = real
+  kind = kind_phys
+  intent = inout
+[pert_scale_plume]
+  standard_name = plume_emissions_scaling_factor
+  long_name = Scaling factor for emissions of plume rising
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[ca_emis_plume]
+  standard_name = fraction_of_cellular_automata_for_plume_rise_emissions
+  long_name = fraction of cellular automata for plume rise emissions
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[biomass_burn_opt_in]
+  standard_name = catchem_biomass_burn_opt
+  long_name = catchem biomass burning option
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[plumerise_flag_in]
+  standard_name = catchem_plumerise_flag
+  long_name = catchem plumerise flag
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[plumerisefire_frq_in]
+  standard_name = catchem_plumerisefire_frq
+  long_name = catchem plumerise fire frequency
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[emis_amp_plume]
+  standard_name = plume_emissions_perturbation_amplitude
+  long_name = multiplier of emissions random perturbation of plume rise emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[ebu]
+  standard_name = chem_buffer_ebu
+  long_name = chemistry buffer of ebu
+  units = various
+  dimensions = (horizontal_loop_extent,vertical_interface_dimension,1,7)
+  type = real
+  kind = kind_phys
+  intent = inout
+[fire_GBBEPx]
+  standard_name = emission_fire_GBBEPx
+  long_name = emission fire GBBEPx
+  units = various
+  dimensions = (horizontal_loop_extent,5)
+  type = real
+  kind = kind_phys
+  intent = in
+[fire_MODIS]
+  standard_name = emission_fire_MODIS
+  long_name = emission fire MODIS
+  units = various
+  dimensions = (horizontal_loop_extent,13)
+  type = real
+  kind = kind_phys
+  intent = in
+# Emissions / Chem input <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<,
+[kemit_in]
+  standard_name = k_emit_dry_dep
+  long_name = catchem k emit dimension for dry deposition
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[pert_scale_anthro]
+  standard_name = anthropogenic_scaling_factor
+  long_name = Scaling factor for emissions of anthropogenic emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[emis_amp_anthro]
+  standard_name = anthropogenic_emissions_perturbation_amplitude
+  long_name = multiplier of emissions random perturbation of anthropogenic emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[emi_in]
+  standard_name = anthropogenic_background_input_cplchp
+  long_name = anthropogenic background input cplchp
+  units = various
+  dimensions = (horizontal_loop_extent,10)
+  type = real
+  kind = kind_phys
+  intent = in
+[dust_in]
+  standard_name = fengsha_dust_input
+  long_name = fengsha dust input
+  units = various
+  dimensions = (horizontal_loop_extent,5)
+  type = real
+  kind = kind_phys
+  intent = in
+[DustFlx]
+  standard_name = instantaneous_dust_emission_flux
+  long_name = instantaneous dust emission flux
+  units = kg m-2 s-1
+  dimensions = (horizontal_loop_extent,5)
+  type = real
+  kind = kind_phys
+# Chemistry Diagnostics <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[drydep]
+  standard_name = instantaneous_dry_deposition
+  long_name = instantaneous dry deposition
+  units = kg m-2 s-1
+  dimensions = (horizontal_loop_extent,number_of_chemical_tracers_for_diagnostics)
+  type = real
+  kind = kind_phys
+  intent = inout
+[wetdpl]
+  standard_name = instantaneous_large_scale_wet_deposition
+  long_name = instantaneous large-scale wet deposition
+  units = kg m-2 s-1
+  dimensions = (horizontal_loop_extent,number_of_chemical_tracers_for_diagnostics)
+  type = real
+  kind = kind_phys
+  intent = inout
+# Cellular Automata <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[ca_sgs_emis]
+  standard_name = flag_for_sgs_cellular_automata_in_chemical_tracer_emissions
+  long_name = switch for sgs ca in chemical tracer emissions
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[ca_sgs]
+  standard_name = flag_for_sgs_cellular_automata
+  long_name = switch for sgs ca
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+# SPPT options <<<<<<<<<<<<<<<<<<<<<<<<<<<<
+[emis_amp_seas]
+  standard_name = sea_spray_emissions_perturbation_amplitude
+  long_name = multiplier of emissions random perturbation of sea salt emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[do_sppt_emis]
+  standard_name = flag_for_stochastic_emissions_perturbations
+  long_name = flag for stochastic emissions perturbations
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[sppt_wts]
+  standard_name = sppt_weights_from_coupled_process
+  long_name = weights for stochastic sppt perturbation
+  units = none
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[ca_sgs_gbbepx_frp]
+  standard_name = GBBEPx_fire_radiative_power_for_stochastic_physics
+  long_name = GBBEPx fire radiative power for stochastic physics
+  units = MW
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = out
+[emis_amp_dust]
+  standard_name = dust_emissions_perturbation_amplitude
+  long_name = multiplier of emissions random perturbation of dust emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[pert_scale_dust]
+  standard_name = dust_emissions_scaling_factor
+  long_name = Scaling factor for emissions of dust emissions
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
diff --git a/drivers/ccpp/ccpp_wrapper.F90 b/drivers/ccpp/ccpp_wrapper.F90
new file mode 100644
index 00000000..478f7023
--- /dev/null
+++ b/drivers/ccpp/ccpp_wrapper.F90
@@ -0,0 +1,792 @@
+!> \file catchem_wrapper_utils.F90
+!! \brief CATCHEM-CCPP interface utilities module
+!!
+!! \details
+!! This is the CCPP-Compliant wrapper for interfacing CATCHEM chemistry model with CCPP
+!! framework. Handles data transformation and management between host model and
+!! CATCHEM chemistry calculations. This module contains the init run and finalize functions and
+!! subroutines that facilitate the data exchange and coordinate transformations
+!! required for CATCHEM integration within CCPP-compliant models.
+!!
+!! \author Barry Baker, NOAA/OAR/ARL
+!!
+!! \date November 2024
+!! \defgroup catchem_ccpp_group CATChem CCPP Interface
+!! \ingroup catchem_ccpp_group
+!!
+!! \note This is part of the CATCHEM-CCPP interface layer that enables
+!!       chemistry calculations within the CCPP framework
+!!!>
+module catchem_wrapper
+
+   use catchem_types, only: catchem_container_type !> CATChem container type
+
+   use physcons, only: g => con_g, pi => con_pi
+   use machine, only: kind_phys
+   use catchem_config
+   use CATChem
+   use catchem_wrapper_utils
+
+   implicit none
+
+   private
+
+   public :: catchem_wrapper_init, catchem_runphase1_wrapper_run
+
+   type(ConfigType) :: Config                          !> CATChem configuration object
+   type(catchem_container_type) :: CATChemStates       !> Container for all CATChem states
+
+!   integer :: im    !> Number of horizontal points
+!   integer :: kme   !> Number of vertical levels
+!   integer :: nsoil !> Number of soil layers
+
+contains
+
+
+
+   !> \section arg_table_catchem_wrapper_init Argument Table
+   !! \htmlinclude catchem_wrapper_init.html
+   !!
+   !! \brief Initialize the CATChem container
+   !! \param[in] catchem_configfile_in Name of the CATChem configuration file
+   !! \param[in] do_catchem_in Flag to enable CATChem calculations
+   !! \param[in] export_catchem_diags_in Flag to export CATChem diagnostics
+   !! \param[in] n_dbg_lines_in Number of debug lines
+   !! \param[in] im Number of horizontal points
+   !! \param[in] kme Number of vertical levels
+   !! \param[in] nsoil Number of soil layers
+   !! \param[out] errmsg Error message
+   !! \param[out] errflg Error flag
+   !!
+   !! \note This subroutine initializes the CATChem container and reads the configuration
+   !!       file to set up the chemistry model
+   !!
+   !! \ingroup catchem_ccpp_group
+   !!!>
+   subroutine catchem_wrapper_init(catchem_configfile_in, do_catchem_in, &
+                                 export_catchem_diags_in, n_dbg_lines_in, &
+                                 im, kme, nsoil, errmsg, errflg)
+
+      use CATChem, only: ConfigType, GridStateType
+      use catchem_wrapper_utils, only: read_catchem_config, allocate_states
+
+      implicit none
+
+      ! Input parameters
+      character(len=*), intent(in) :: catchem_configfile_in
+      logical,          intent(in) :: do_catchem_in
+      logical,          intent(in) :: export_catchem_diags_in
+      integer,          intent(in) :: n_dbg_lines_in
+      integer,          intent(in) :: im
+      integer,          intent(in) :: kme
+      integer,          intent(in) :: nsoil
+
+      ! Output parameters
+      character(len=*), intent(out) :: errmsg
+      integer,          intent(out) :: errflg
+
+      ! Local variables
+      integer :: i
+
+      errmsg = ''
+      errflg = 0
+
+      if (.not. do_catchem_in) return
+
+      ! Set global parameters
+      CATCHem_ConfigFile = catchem_config_opt
+      do_catchem = do_catchem_in
+      export_catchem_diags = export_catchem_diags_in
+      n_dbg_lines = n_dbg_lines_in
+
+      ! Initialize CATChem container
+      call CATChemStates%init(im)
+
+      ! Read configuration
+      call read_catchem_config(Config, CATChemStates, CATCHem_ConfigFile, errflg, errmsg)
+      if (errflg /= 0) return
+
+
+
+      ! Allocate states for each horizontal point
+      do i = 1, im
+         call allocate_states(Config, MetState(i), ChemState(i), DiagState(i), &
+                              EmisState(i), kme, nsoil, errflg, errmsg)
+         if (errflg /= 0) return
+      end do
+
+   end subroutine catchem_wrapper_init
+
+  !> \brief Brief description of the subroutine
+  !!
+  !! \section arg_table_catchem_gocart_wrapper_finalize Argument Table
+  !!
+  subroutine catchem_wrapper_finalize(im, kme, nsoil, errmsg, errflg)
+
+    use catchem_wrapper_utils, only: deallocate_states
+
+   implicit none
+
+   ! Input parameters
+   integer, intent(in) :: im
+   integer, intent(in) :: kme
+   integer, intent(in) :: nsoil
+
+   ! Output parameters
+   character(len=*), intent(out) :: errmsg
+   integer, intent(out) :: errflg
+
+   ! Local variables
+   integer :: i !> Horizontal index
+
+   do i = 1, im
+      call deallocate_states(MetState(i), ChemState(i), DiagState(i), EmisState(i), kme, nsoil, errflg, errmsg)
+      if (errflg /= 0) return
+   end do
+
+   deallocate(CATChemStates)
+
+  end subroutine catchem_wrapper_finalize
+
+  !>
+  !! This is the Configurable ATmospheric Chemistry (CATChem)
+  !! This is the CATChem gocart wrapper Module
+  !! \section arg_table_catchem_runphase1_wrapper_run Argument Table
+  !! \htmlinclude catchem_runphase1_wrapper_run.html
+  !!
+  !>\section catchem_phase1_group CATChem Scheme General Algorithm
+  !> @{
+  subroutine catchem_wrapper_run(im, kte, kme, ktau, dt, jdate, julian, idat, &
+     garea, rlat, rlon, xcosz, &
+     ntrac, ntc, ntr, tile_num, mpicomm, mpirank, mpiroot, &
+     land, vegtype_dom, nlcat, vegtype_frac, oro, &
+     nsoil, smc, soiltyp, tslb, snow, pgr, pb2d, hpbl_thetav, kpbl, &
+     sncovr1_lnd, albedo, &
+     u10m, v10m, ustar, tskin, t2m, dpt2m, hf2d, lf2d,  &
+     znt, dswsfc, recmol, rain_cpl, rainc_cpl, &
+     pr3d, ph3d, phl3d, prl3d, tk3d, us3d, vs3d, w, prslp, q3d, &
+     epsm1,   &
+     qrn, qsnw, &
+     cldf, &
+     qg0, qgrs, nwfa, nifa, &
+     gq0, qgrs, errmsg, errflg)
+
+     implicit none
+
+     ! ARGUMENTS
+     !----------
+
+     ! Grid information
+     integer, intent(in) :: im       ! horizontal loop extent
+     integer, intent(in) :: kte      ! vertical layer dimension
+     integer, intent(in) :: kme      ! vertical interface dimension = number of vertical layers + 1
+     integer, intent(in) :: ktau     ! current forecast iteration
+     integer, intent(in) :: tile_num ! index of cube sphere tile
+     real(kind_phys), dimension(im), intent(in) :: garea !> grid area (m^2) of each grid cell
+     real(kind_phys), dimension(im), intent(in) :: rlat  !> latitude
+     real(kind_phys), dimension(im), intent(in) :: rlon, !> longitude
+     real(kind_phys), dimension(im), intent(in) :: xcosz !> cosine of solar zenith angle
+
+
+     ! Time information
+     integer, intent(in) :: idat(8)  ! initialization date and time (in iso order)
+     integer, intent(in) :: jdate    !> julian date
+     real(kind_phys), intent(in) :: dt      ! physics timestep (s)
+     real(kind_phys), intent(in) :: julian  ! forecast julian day (days)
+
+     ! MPI information
+     type(MPI_comm), intent(in) :: mpicomm
+     integer, intent(in) :: mpirank
+     integer, intent(in) :: mpiroot
+
+     ! tracer information
+     integer, intent(in) :: ntrac    ! total number of tracers
+     integer, intent(in) :: ntc      ! number of chemical tracers
+     integer, intent(in) :: ntr      ! number of aerosol tracers
+     real(kind_phys), dimension(im, kte, ntrac), intent(inout) :: gq0
+     real(kind_phys), dimension(im, kte, ntrac), intent(inout) :: qgrs
+
+     ! land surface information
+     integer, dimension(im), intent(in)                :: land          !> sea land ice mask (sea = 0, land = 1, ice = 2)
+     integer, dimension(im), intent(in)                :: soiltyp       !> soil type
+     integer, dimension(im), intent(in)                :: vegtype_dom   !> vegetation type
+     integer,                intent(in)                :: nlcat         !> number of land surface categories
+     real(kind_phys), dimension(im, nlcat), intent(in) :: vegtype_frac  !> fraction of each land surface category
+     real(kind_phys), dimension(im), intent(in)        :: oro           !> height above mean sea level (m)
+     real(kind_phys), dimension(im), intent(in)        :: nsoil         !> number of soil layers
+     real(kind_phys), dimension(im, nsoil), intent(in) :: smc           !> volumetric fraction of soil moisture for lsm
+     real(kind_phys), dimension(im,nsoil), intent(in)  :: tslb          !> soil temperature (K
+     real(kind_phys), dimension(im), intent(in)        :: snow          !> water equivalent snow depth (mm)
+     real(kind_phys), dimension(im), intent(in)        :: sncovr1_lnd   !> fractional snow cover land
+     real(kind_phys), dimension(im), intent(in)        :: pgr           !> pressure at the surface (Pa)
+     real(kind_phys), dimension(im), intent(in)        :: pb2d          !> PBL Thickness determined by the PBL scheme (m)
+     real(kind_phys), dimension(im), intent(in)        :: kpbl          !> PBL level
+     real(kind_phys), dimension(im), intent(in)        :: hpbl_thetav   !> PBL Height based on modified parcel method (m)
+     real(kind_phys), dimension(im), intent(in)        :: u10m          !> 10 m wind speed (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: v10m          !> 10 m wind speed (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: ustar         !> friction velocity (m/s)
+     real(kind_phys), dimension(im), intent(in)        :: tskin         !> skin temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: t2m           !> 2 m temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: dpt2m         !> 2 m dew point temperature (K)
+     real(kind_phys), dimension(im), intent(in)        :: hf2d          !> Sensible heat flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: lf2d          !> Latent heat flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: znt          !> surface roughness length in (cm)
+     real(kind_phys), dimension(im), intent(in)        :: dswsfc        !> downward short wave flux (W m-2)
+     real(kind_phys), dimension(im), intent(in)        :: recmol        !> one over obukhov length (m-1)
+     real(kind_phys), dimension(im), intent(in)        :: albedo        !> surface albedo
+     real(kind_phys), dimension(im), intent(in)        :: prslp         !> sea level pressure (Pa)
+
+     ! 3d state information
+     real(kind_phys), dimension(im, kte), intent(in) :: pr3d            !> air pressure at model layer interfaces (Pa)
+     real(kind_phys), dimension(im, kte), intent(in) :: prl3d           !> pressure at the model level (Pa)
+     real(kind_phys), dimension(im, kte), intent(in) :: ph3d            !> geopotential at the model level (m2 s-2)
+     real(kind_phys), dimension(im, kte), intent(in) :: phl3d           !> geopotential at the model layer interfaces (m2 s-2)
+     real(kind_phys), dimension(im, kte), intent(in) :: tk3d
+     real(kind_phys), dimension(im, kte), intent(in) :: us3d
+     real(kind_phys), dimension(im, kte), intent(in) :: vs3d
+     real(kind_phys), dimension(im, kte), intent(in) :: q3d
+
+     ! precipitation information
+     real(kind_phys), dimension(im), intent(in)        :: rain_cpl      !> total rain at this time step (m)
+     real(kind_phys), dimension(im), intent(in)        :: rainc_cpl     !> convective rain at this time step (m)
+     real(kind_phys), dimension(im), intent(in)        :: cldf          !> total cloud fraction
+
+     real(kind_phys), dimension(im, kte, 3), intent(in) :: emi2_in ! 3 should be temporary... need to replace with either an input namelist option or have it be a pointer with variable dimension
+
+     ! Output
+     character(len=*), intent(out) :: errmsg
+     integer, intent(out) :: errflg
+     real(kind_phys), dimension(im, kte, 3), intent(out) :: emi2_out ! 3 should be temporary... need to replace with either an input namelist option or have it be a pointer with variable dimension
+     ! Local
+     integer :: mpiid
+     integer, parameter :: ids = 1, jds = 1, jde = 1, kds = 1
+     integer, parameter :: ims = 1, jms = 1, jme = 1, kms = 1
+     integer, parameter :: its = 1, jts = 1, jte = 1, kts = 1
+
+     real(kind_phys), dimension(1:im, 1:kme, jms:jme) :: rri, t_phy, &
+        p_phy, dz8w, p8w, t8w, rho_phy
+
+     real(kind_phys), dimension(ims:im, jms:jme) :: xlat, xlong, dxy
+
+     !>- chemistry variables
+     real(kind_phys), dimension(ims:im, kms:kme, jms:jme, 1:num_moist)  :: moist
+     real(kind_phys), dimension(ims:im, kms:kme, jms:jme, 1:num_chem)  :: chem
+
+     integer :: ide, ime, ite, kde, julday
+
+     !   integer, parameter :: SEAS_OPT_DEFAULT = 1
+     !   integer, parameter :: chem_in_opt = 0  ! 0 for coldstart, 1 for restart
+     logical, parameter :: readrestart = .false.
+     integer, parameter :: nvl_gocart = 64  ! number of input levels from gocart file
+
+     real(kind_phys), dimension(ims:im, kms:kme, jms:jme) :: pm10, pm2_5_dry, pm2_5_dry_ec
+
+     !>- chemical background variables
+     real(kind_phys), dimension(ims:im, kms:kme, jms:jme) :: backg_oh, backg_h2o2, backg_no3
+
+     real(kind_phys), dimension(ims:im, kms:kme, jms:jme) :: oh_t, h2o2_t, no3_t
+     real(kind_phys), dimension(ims:im, jms:jme) :: ttday, tcosz
+
+     real(kind_phys) :: dtstep, gmt
+     real(kind_phys), dimension(1:num_chem) :: ppm2ugkg
+
+     ! -- output tracers
+     real(kind_phys), dimension(ims:im, jms:jme, 1:kme) :: p10, pm25!, ebu_oc
+     real(kind_phys), dimension(ims:im, jms:jme, 1:kme) :: oh_bg, h2o2_bg, no3_bg
+
+     !>-- local variables
+     logical :: call_gocart
+     integer :: i, j, jp, k, kp, n
+
+     errmsg = ''
+     errflg = 0
+
+     chem_opt = chem_opt_in
+
+     gmt = real(idat(5))
+     julday = real(julian)
+
+     ! -- set domain
+     ide = im
+     ime = im
+     ite = im
+     kde = kte
+
+     ! -- volume to mass fraction conversion table (ppm -> ug/kg)
+     ppm2ugkg = 1._kind_phys
+     ppm2ugkg(p_sulf) = 1.e+03_kind_phys*mw_so4_aer/mwdry
+
+     ! -- set control flags
+     call_gocart = (mod(ktau, call_chemistry) == 0) .or. (ktau == 1)
+
+     ! -- compute accumulated large-scale and convective rainfall since last call
+     if (ktau > 1) then
+        dtstep = call_chemistry*dt
+     else
+        dtstep = dt
+     end if
+
+     !!!
+
+     !>- get ready for chemistry run
+     call catchem_gocart_prep( &
+        readrestart, chem_in_opt, ktau, dtstep, xcosz, &
+        garea, rlat, rlon, &
+        pr3d, ph3d, tk3d, prl3d, spechum, &
+        emi2_in, &
+        xlat, xlong, dxy, &
+        rri, t_phy, p_phy, rho_phy, dz8w, p8w, &
+        t8w, &
+        ntso2, ntsulf, ntDMS, ntmsa, ntpp25, &
+        ntbc1, ntbc2, ntoc1, ntoc2, ntpp10, &
+        ntrac, gq0, &
+        num_chem, num_moist, &
+        call_gocart, nvl_gocart, &
+        ttday, tcosz, gmt, julday, &
+        backg_oh, backg_h2o2, backg_no3, &
+        ppm2ugkg, &
+        moist, chem, &
+        ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte)
+     !write (*,*) 'hli test2 ktau',call_gocart
+
+     if (call_gocart) then
+        call gocart_chem_driver(ktau, dt, dtstep, gmt, julday, xcosz, &
+           t_phy, moist, chem, rho_phy, dz8w, p8w, backg_oh, oh_t, &
+           backg_h2o2, h2o2_t, backg_no3, no3_t, &
+           dxy, g, xlat, xlong, ttday, tcosz, &
+           chem_opt, num_chem, num_moist, &
+           ids, ide, jds, jde, kds, kde, &
+           ims, ime, jms, jme, kms, kme, &
+           its, ite, jts, jte, kts, kte)
+        call gocart_aerosols_driver(ktau, dtstep, t_phy, moist, &
+           chem, rho_phy, dz8w, p8w, dxy, g, &
+           chem_opt, num_chem, num_moist, &
+           ids, ide, jds, jde, kds, kde, &
+           ims, ime, jms, jme, kms, kme, &
+           its, ite, jts, jte, kts, kte)
+     end if
+
+     call sum_pm_gocart( &
+        rri, chem, pm2_5_dry, pm2_5_dry_ec, pm10, &
+        num_chem, chem_opt, &
+        ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte)
+
+     ! -- pm25 and pm10 for output , not for tracer options
+     do j = jts, jte
+        do k = kts, kte
+           do i = its, ite
+              pm25(i, j, k) = pm2_5_dry(i, k, j)
+              p10(i, j, k) = pm10(i, k, j)
+              !ebu_oc(i,j,k) = ebu      (i,k,j,p_ebu_oc)
+           end do
+        end do
+     end do
+
+     if (call_gocart) then
+        do j = jts, jte
+           do k = kts, kte
+              do i = its, ite
+                 oh_bg(i, j, k) = max(0., oh_t(i, k, j))
+                 h2o2_bg(i, j, k) = max(0., h2o2_t(i, k, j))
+                 no3_bg(i, j, k) = max(0., no3_t(i, k, j))
+              end do
+           end do
+        end do
+     end if
+
+     ! -- put chem stuff back into tracer array
+     do k = kts, kte
+        do i = its, ite
+           gq0(i, k, ntso2) = ppm2ugkg(p_so2)*max(epsilc, chem(i, k, 1, p_so2))
+           gq0(i, k, ntsulf) = ppm2ugkg(p_sulf)*max(epsilc, chem(i, k, 1, p_sulf))
+           gq0(i, k, ntdms) = ppm2ugkg(p_dms)*max(epsilc, chem(i, k, 1, p_dms))
+           gq0(i, k, ntmsa) = ppm2ugkg(p_msa)*max(epsilc, chem(i, k, 1, p_msa))
+           gq0(i, k, ntpp25) = ppm2ugkg(p_p25)*max(epsilc, chem(i, k, 1, p_p25))
+           gq0(i, k, ntbc1) = ppm2ugkg(p_bc1)*max(epsilc, chem(i, k, 1, p_bc1))
+           gq0(i, k, ntbc2) = ppm2ugkg(p_bc2)*max(epsilc, chem(i, k, 1, p_bc2))
+           gq0(i, k, ntoc1) = ppm2ugkg(p_oc1)*max(epsilc, chem(i, k, 1, p_oc1))
+           gq0(i, k, ntoc2) = ppm2ugkg(p_oc2)*max(epsilc, chem(i, k, 1, p_oc2))
+           gq0(i, k, ntpp10) = ppm2ugkg(p_p10)*max(epsilc, chem(i, k, 1, p_p10))
+        end do
+     end do
+
+     do k = kts, kte
+        do i = its, ite
+           qgrs(i, k, ntso2) = gq0(i, k, ntso2)
+           qgrs(i, k, ntsulf) = gq0(i, k, ntsulf)
+           qgrs(i, k, ntdms) = gq0(i, k, ntdms)
+           qgrs(i, k, ntmsa) = gq0(i, k, ntmsa)
+           qgrs(i, k, ntpp25) = gq0(i, k, ntpp25)
+           qgrs(i, k, ntbc1) = gq0(i, k, ntbc1)
+           qgrs(i, k, ntbc2) = gq0(i, k, ntbc2)
+           qgrs(i, k, ntoc1) = gq0(i, k, ntoc1)
+           qgrs(i, k, ntoc2) = gq0(i, k, ntoc2)
+           qgrs(i, k, ntpp10) = gq0(i, k, ntpp10)
+        end do
+     end do
+
+     !
+  end subroutine catchem_gocart_wrapper_run
+  !> @}
+
+  subroutine catchem_gocart_prep( &
+     readrestart, chem_in_opt, ktau, dtstep, xcosz, &
+     garea, rlat, rlon, &
+     pr3d, ph3d, tk3d, prl3d, spechum, &
+     emi2_in, &
+     xlat, xlong, dxy, &
+     rri, t_phy, p_phy, rho_phy, dz8w, p8w, &
+     t8w, &
+     ntso2, ntsulf, ntDMS, ntmsa, ntpp25, &
+     ntbc1, ntbc2, ntoc1, ntoc2, ntpp10, &
+     ntrac, gq0, &
+     num_chem, num_moist, &
+     call_gocart, nvl_gocart, &
+     ttday, tcosz, gmt, julday, &
+     backg_oh, backg_h2o2, backg_no3, &
+     ppm2ugkg, &
+     moist, chem, &
+     ids, ide, jds, jde, kds, kde, &
+     ims, ime, jms, jme, kms, kme, &
+     its, ite, jts, jte, kts, kte)
+
+     !Chem input configuration
+     logical, intent(in) :: readrestart
+     integer, intent(in) :: chem_in_opt, ktau, julday
+     real(kind=kind_phys), intent(in) :: dtstep, gmt
+
+     !FV3 input variables
+     integer, intent(in) :: ntrac
+     integer, intent(in) :: ntso2, ntpp25, ntbc1, ntoc1, ntpp10
+     integer, intent(in) :: ntsulf, ntbc2, ntoc2, ntDMS, ntmsa
+     real(kind=kind_phys), dimension(ims:ime), intent(in) :: garea, rlat, rlon, xcosz
+     real(kind=kind_phys), dimension(ims:ime, 64, 3), intent(in) :: emi2_in
+     real(kind=kind_phys), dimension(ims:ime, kms:kme), intent(in) :: pr3d, ph3d
+     real(kind=kind_phys), dimension(ims:ime, kts:kte), intent(in) :: &
+        tk3d, prl3d, spechum
+     real(kind=kind_phys), dimension(ims:ime, kts:kte, ntrac), intent(in) :: gq0
+
+     !GSD Chem variables
+     integer, intent(in) ::  num_chem, num_moist, &
+        nvl_gocart
+     logical, intent(in) ::  call_gocart
+     integer, intent(in) ::  ids, ide, jds, jde, kds, kde, &
+        ims, ime, jms, jme, kms, kme, &
+        its, ite, jts, jte, kts, kte
+
+     real(kind_phys), dimension(num_chem), intent(in) :: ppm2ugkg
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: &
+        backg_oh, backg_h2o2, backg_no3
+
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: &
+        rri, t_phy, p_phy, rho_phy, dz8w, p8w, t8w
+     real(kind_phys), dimension(ims:ime, jms:jme), intent(out) :: &
+        xlat, xlong, dxy, &
+        ttday, tcosz
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_moist), intent(out) :: moist
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_chem), intent(out) :: chem
+
+     real(kind_phys), dimension(ims:ime, kms:kme, jms:jme) :: z_at_w
+     real(kind_phys), dimension(nvl_gocart) :: p_gocart
+
+     ! -- local variables
+     real(kind_phys), dimension(ims:ime, jms:jme, nvl_gocart) :: oh_backgd, h2o2_backgd, no3_backgd
+     real(kind_phys) ::  pu, pl, aln, pwant
+     real(kind_phys) ::  xhour, xmin, gmtp, xlonn, xtime, real_time
+     real(kind_phys), DIMENSION(1, 1) :: sza, cosszax
+     integer i, ip, j, jp, k, kp, kk, kkp, nv, jmax, jmaxi, l, ll, n, ndystep, ixhour
+
+     p_gocart = (/1000., 992.5, 985., 977.5, 970., 955., 940., 925., 910., &
+        895., 880., 865., 850., 825., 800., 775., 750., 712.5, 675., 637.5, 600., &
+        562.5, 525., 487.5, 450., 412.5, 375., 337.5, 288.08, 244.88, 208.15, 176.93, &
+        150.39, 127.84, 108.66, 92.37, 78.51, 66.6, 56.39, 47.64, 40.18, 33.81, 28.37, &
+        23.73, 19.79, 16.46, 13.64, 11.28, 9.29, 7.62, 6.22, 5.05, 4.08, 3.28, 2.62, &
+        2.08, 1.65, 1.3, 1.02, 0.8, 0.62, 0.48, 0.37, 0.28/)
+
+     ! -- initialize output arrays
+     backg_oh = 0._kind_phys
+     backg_h2o2 = 0._kind_phys
+     backg_no3 = 0._kind_phys
+     rri = 0._kind_phys
+     t_phy = 0._kind_phys
+     p_phy = 0._kind_phys
+     rho_phy = 0._kind_phys
+     dz8w = 0._kind_phys
+     p8w = 0._kind_phys
+     t8w = 0._kind_phys
+     xlat = 0._kind_phys
+     xlong = 0._kind_phys
+     dxy = 0._kind_phys
+     ttday = 0._kind_phys
+     tcosz = 0._kind_phys
+     moist = 0._kind_phys
+     chem = 0._kind_phys
+     z_at_w = 0._kind_phys
+
+     do i = its, ite
+        dxy(i, 1) = garea(i)
+        xlat(i, 1) = rlat(i)*180./pi
+        xlong(i, 1) = rlon(i)*180./pi
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do i = its, ite
+           ip = i - its + 1
+           z_at_w(i, kts, j) = max(0., ph3d(ip, 1)/g)
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte
+           kp = k - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              dz8w(i, k, j) = abs(ph3d(ip, kp + 1) - ph3d(ip, kp))/g
+              z_at_w(i, k + 1, j) = z_at_w(i, k, j) + dz8w(i, k, j)
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte + 1
+           kp = k - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              p8w(i, k, j) = pr3d(ip, kp)
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        jp = j - jts + 1
+        do k = kts, kte + 1
+           kk = min(k, kte)
+           kkp = kk - kts + 1
+           do i = its, ite
+              ip = i - its + 1
+              dz8w(i, k, j) = z_at_w(i, kk + 1, j) - z_at_w(i, kk, j)
+              t_phy(i, k, j) = tk3d(ip, kkp)
+              p_phy(i, k, j) = prl3d(ip, kkp)
+              rho_phy(i, k, j) = p_phy(i, k, j)/(287.04*t_phy(i, k, j)*(1.+.608*spechum(ip, kkp)))
+              rri(i, k, j) = 1./rho_phy(i, k, j)
+              moist(i, k, j, :) = 0.
+              moist(i, k, j, 1) = gq0(ip, kkp, p_atm_shum)
+              if (t_phy(i, k, j) > 265.) then
+                 moist(i, k, j, 2) = gq0(ip, kkp, p_atm_cldq)
+                 moist(i, k, j, 3) = 0.
+                 if (moist(i, k, j, 2) < 1.e-8) moist(i, k, j, 2) = 0.
+              else
+                 moist(i, k, j, 2) = 0.
+                 moist(i, k, j, 3) = gq0(ip, kkp, p_atm_cldq)
+                 if (moist(i, k, j, 3) < 1.e-8) moist(i, k, j, 3) = 0.
+              end if
+              !--
+           end do
+        end do
+     end do
+
+     do j = jts, jte
+        do k = 2, kte
+           do i = its, ite
+              t8w(i, k, j) = .5*(t_phy(i, k, j) + t_phy(i, k - 1, j))
+           end do
+        end do
+     end do
+
+     ! -- only used in phtolysis....
+     do j = jts, jte
+        do i = its, ite
+           t8w(i, 1, j) = t_phy(i, 1, j)
+           t8w(i, kte + 1, j) = t_phy(i, kte, j)
+        end do
+     end do
+
+     do k = kms, kte
+        do i = ims, ime
+           chem(i, k, jts, p_so2) = max(epsilc, gq0(i, k, ntso2)/ppm2ugkg(p_so2))
+           chem(i, k, jts, p_sulf) = max(epsilc, gq0(i, k, ntsulf)/ppm2ugkg(p_sulf))
+           chem(i, k, jts, p_dms) = max(epsilc, gq0(i, k, ntdms)/ppm2ugkg(p_dms))
+           chem(i, k, jts, p_msa) = max(epsilc, gq0(i, k, ntmsa)/ppm2ugkg(p_msa))
+           chem(i, k, jts, p_p25) = max(epsilc, gq0(i, k, ntpp25)/ppm2ugkg(p_p25))
+           chem(i, k, jts, p_bc1) = max(epsilc, gq0(i, k, ntbc1)/ppm2ugkg(p_bc1))
+           chem(i, k, jts, p_bc2) = max(epsilc, gq0(i, k, ntbc2)/ppm2ugkg(p_bc2))
+           chem(i, k, jts, p_oc1) = max(epsilc, gq0(i, k, ntoc1)/ppm2ugkg(p_oc1))
+           chem(i, k, jts, p_oc2) = max(epsilc, gq0(i, k, ntoc2)/ppm2ugkg(p_oc2))
+           chem(i, k, jts, p_p10) = max(epsilc, gq0(i, k, ntpp10)/ppm2ugkg(p_p10))
+        end do
+     end do
+
+     if (.NOT. readrestart) then
+        if (chem_in_opt == 0) then
+           if (ktau .le. 1) then
+              !           if(chem_opt > 0 ) then
+              do j = jts, jte
+                 jp = j - jts + 1
+                 do k = kts, kte
+                    do i = its, ite
+                       ip = i - its + 1
+                       if (chem_opt == CHEM_OPT_GOCART) then
+                          do n = 1, num_chem
+                             chem(i, k, j, n) = 1.e-30
+                          end do
+                       end if  ! chem_opt==300
+                       if ((chem_opt > CHEM_OPT_GOCART) .and. (chem_opt < CHEM_OPT_MAX)) then
+                          chem(i, k, j, p_so2) = 5.e-10
+                          chem(i, k, j, p_sulf) = 3.e-10
+                          chem(i, k, j, p_msa) = 1.e-10
+                          chem(i, k, j, p_dms) = 1.e-10
+                       end if !chem_opt >= 300 .and. chem_opt <  500
+
+                       !                if ((chem_opt == CHEM_OPT_GOCART_RACM) .or. (chem_opt == CHEM_OPT_RACM_SOA_VBS)) then  !added o3 background !lzhang
+                       !                  kk=min(k,kte)
+                       !                  kkp = kk - kts + 1
+                       !                  ! -- add initial constant into O3,CH4 and CO ect.
+                       !                  chem(i,k,j,p_o3)=epsilc
+                       !                  ! -- this section needs to be revisited before enabling the
+                       !                  ! corresponding chem_opt options
+                       !                  ! maxth=min(400.,th_pvsrf(i,j))
+                       !                  ! if (tr3d(ip,jp,kkp,p_atm_ptem) > maxth) then
+                       !                  !   chem(i,k,j,p_o3)=(airmw/48.)*tr3d(ip,jp,kkp,p_atm_o3mr)*1e6
+                       !                  !   !convert kg/kg to ppm
+                       !                  ! else
+                       !                  !   chem(i,k,j,p_o3)=0.03 !ppm
+                       !                  ! endif
+                       !                  chem(i,k,j,p_ch4)=1.85 !ppm
+                       !                  chem(i,k,j,p_co)=0.06 !ppm
+                       !                  chem(i,k,j,p_co2)=380.
+                       !                  chem(i,k,j,p_ete)=epsilc
+                       !                  chem(i,k,j,p_udd)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_hket)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_api)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_lim)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_dien)=chem(i,k,j,p_ete)
+                       !                  chem(i,k,j,p_macr)=chem(i,k,j,p_ete)
+                       !                endif !( (chem_opt == 301.or.chem_opt==108))
+                    end do
+                 end do
+              end do
+           end if !(ktau<=1)
+
+        else !(chem_in_opt == 0 )
+
+           if ((ktau <= 1) .and. ((chem_opt == CHEM_OPT_GOCART_RACM) .or. (chem_opt == CHEM_OPT_RACM_SOA_VBS))) then  !added GFS o3 background above 380K!lzhang
+              do j = jts, jte
+                 jp = j - jts + 1
+                 do k = kts, kte + 1
+                    kk = min(k, kte)
+                    kkp = kk - kts + 1
+                    do i = its, ite
+                       ip = i - its + 1
+                       ! -- this section needs to be revisited before enabling the
+                       ! corresponding chem_opt options
+                       ! maxth=min(400.,th_pvsrf(i,j))
+                       ! if (tr3d(ip,jp,kkp,p_atm_ptem) >= maxth) then
+                       !   chem(i,k,j,p_o3)=(airmw/48.)*tr3d(ip,jp,kkp,p_atm_o3mr)*1e6 !convert kg/kg to ppm
+                       ! endif !380K
+                    end do
+                 end do
+              end do
+           end if ! chem_opt == 301.or.chem_opt==108
+
+        end if !(chem_in_opt == 1 )
+     end if ! readrestart
+
+     !-- assgin read in 3D background chemical species
+     do i = its, ite
+        do k = 1, nvl_gocart
+           h2o2_backgd(i, 1, k) = emi2_in(i, k, 1)
+           no3_backgd(i, 1, k) = emi2_in(i, k, 2)
+           oh_backgd(i, 1, k) = emi2_in(i, k, 3)
+        end do
+     end do
+
+     !
+     ! -- gocart background fields only if gocart is called
+     if (call_gocart .and. (chem_opt == CHEM_OPT_GOCART)) then
+        do j = jts, jte
+           do i = its, ite
+              do k = kts, kte
+                 do ll = 2, nvl_gocart
+                    l = ll
+                    if (p_gocart(l) < .01*p_phy(i, k, j)) exit
+                 end do
+                 pu = alog(p_gocart(l))
+                 pl = alog(p_gocart(l - 1))
+                 pwant = alog(.01*p_phy(i, k, j))
+                 if (pwant > pl) then
+                    backg_oh(i, k, j) = oh_backgd(i, j, l)
+                    backg_h2o2(i, k, j) = h2o2_backgd(i, j, l)
+                    backg_no3(i, k, j) = no3_backgd(i, j, l)
+                 else
+                    aln = (oh_backgd(i, j, l)*(pwant - pl) + &
+                       oh_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_oh(i, k, j) = aln
+                    aln = (h2o2_backgd(i, j, l)*(pwant - pl) + &
+                       h2o2_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_h2o2(i, k, j) = aln
+                    aln = (no3_backgd(i, j, l)*(pwant - pl) + &
+                       no3_backgd(i, j, l - 1)*(pu - pwant))/(pu - pl)
+                    backg_no3(i, k, j) = aln
+                 end if
+              end do
+           end do
+        end do
+     end if   ! end gocart stuff
+     !endif !restart
+
+     if ((chem_opt == CHEM_OPT_RACM_SOA_VBS) .or. (chem_opt >= CHEM_OPT_GOCART .and. chem_opt < CHEM_OPT_MAX)) then
+        !ndystep=86400/ifix(dtstepc)
+        ndystep = 86400/ifix(dtstep)
+        do j = jts, jte
+           do i = its, ite
+              tcosz(i, j) = 0.
+              ttday(i, j) = 0.
+              xlonn = xlong(i, j)
+              do n = 1, ndystep
+                 xtime = n*dtstep/60.
+                 ixhour = ifix(gmt + .01) + ifix(xtime/60.)
+                 xhour = float(ixhour)
+                 xmin = 60.*gmt + (xtime - xhour*60.)
+                 gmtp = mod(xhour, 24.)
+                 gmtp = gmtp + xmin/60.
+                 CALL szangle(1, 1, julday, gmtp, sza, cosszax, xlonn, rlat(i))
+                 TCOSZ(i, j) = TCOSZ(I, J) + cosszax(1, 1)
+                 if (cosszax(1, 1) > 0.) ttday(i, j) = ttday(i, j) + dtstep
+              end do
+           end do
+        end do
+     end if !chem_opt >= 300 .and. chem_opt <  500
+
+  end subroutine catchem_gocart_prep
+
+  !> \brief Convert CCPP to CATChem
+  !! \ingroup catchem_wrapper_group
+  subroutine ccpp_to_catchem(chem_in_opt, ktau, dtstep, xcosz, &
+    garea, rlat, rlon, &
+    pr3d, ph3d, tk3d, prl3d, spechum, &
+    emi2_in, &
+    xlat, xlong, dxy, &
+    rri, t_phy, p_phy, rho_phy, dz8w, p8w, &
+    t8w, &
+    ntso2, ntsulf, ntDMS, ntmsa, ntpp25, &
+    ntbc1, ntbc2, ntoc1, ntoc2, ntpp10, &
+    ntrac, gq0, &
+    num_chem, num_moist, &
+    call_gocart, nvl_gocart, &
+    ttday, tcosz, gmt, julday, &
+    backg_oh, backg_h2o2, backg_no3, &
+    ppm2ugkg, &
+    moist, chem, &
+    ids, ide, jds, jde, kds, kde, &
+    ims, ime, jms, jme, kms, kme, &
+    its, ite, jts, jte, kts, kte)
+
+  !> @}
+end module catchem_gocart_wrapper
diff --git a/drivers/ccpp/catchem_anthropogenic_wrapper.F90 b/drivers/ccpp/old/catchem_anthropogenic_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_anthropogenic_wrapper.F90
rename to drivers/ccpp/old/catchem_anthropogenic_wrapper.F90
diff --git a/drivers/ccpp/catchem_anthropogenic_wrapper.meta b/drivers/ccpp/old/catchem_anthropogenic_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_anthropogenic_wrapper.meta
rename to drivers/ccpp/old/catchem_anthropogenic_wrapper.meta
diff --git a/drivers/ccpp/catchem_diag_wrapper.F90 b/drivers/ccpp/old/catchem_diag_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_diag_wrapper.F90
rename to drivers/ccpp/old/catchem_diag_wrapper.F90
diff --git a/drivers/ccpp/catchem_diag_wrapper.meta b/drivers/ccpp/old/catchem_diag_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_diag_wrapper.meta
rename to drivers/ccpp/old/catchem_diag_wrapper.meta
diff --git a/drivers/ccpp/catchem_dmsemis_wrapper.F90 b/drivers/ccpp/old/catchem_dmsemis_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_dmsemis_wrapper.F90
rename to drivers/ccpp/old/catchem_dmsemis_wrapper.F90
diff --git a/drivers/ccpp/catchem_dmsemis_wrapper.meta b/drivers/ccpp/old/catchem_dmsemis_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_dmsemis_wrapper.meta
rename to drivers/ccpp/old/catchem_dmsemis_wrapper.meta
diff --git a/drivers/ccpp/catchem_drydep_wrapper.F90 b/drivers/ccpp/old/catchem_drydep_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_drydep_wrapper.F90
rename to drivers/ccpp/old/catchem_drydep_wrapper.F90
diff --git a/drivers/ccpp/catchem_drydep_wrapper.meta b/drivers/ccpp/old/catchem_drydep_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_drydep_wrapper.meta
rename to drivers/ccpp/old/catchem_drydep_wrapper.meta
diff --git a/drivers/ccpp/catchem_dust_wrapper.F90 b/drivers/ccpp/old/catchem_dust_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_dust_wrapper.F90
rename to drivers/ccpp/old/catchem_dust_wrapper.F90
diff --git a/drivers/ccpp/catchem_dust_wrapper.meta b/drivers/ccpp/old/catchem_dust_wrapper.meta
similarity index 93%
rename from drivers/ccpp/catchem_dust_wrapper.meta
rename to drivers/ccpp/old/catchem_dust_wrapper.meta
index fd368592..c05e466c 100644
--- a/drivers/ccpp/catchem_dust_wrapper.meta
+++ b/drivers/ccpp/old/catchem_dust_wrapper.meta
@@ -68,6 +68,13 @@
   dimensions = (horizontal_loop_extent)
   type = integer
   intent = in
+[nsoil]
+  standard_name = vertical_dimension_of_soil_internal_to_land_surface_scheme
+  long_name = number of soil layers internal to land surface model
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [u10m]
   standard_name = x_wind_at_10m
   long_name = 10 meter u wind speed
@@ -92,6 +99,38 @@
   type = real
   kind = kind_phys
   intent = in
+[pblh]
+  standard_name = atmosphere_boundary_layer_thickness
+  long_name = PBL thickness
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[tskin]
+  standard_name = surface_skin_temperature
+  long_name = surface skin temperature
+  units = K
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[t2m]
+  standard_name = air_temperature_at_2m
+  long_name = 2 meter temperature
+  units = K
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[dpt2m]
+  standard_name = dewpoint_temperature_at_2m
+  long_name = 2 meter dewpoint temperature
+  units = K
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
 [rlat]
   standard_name = latitude
   long_name = latitude
@@ -108,14 +147,7 @@
   type = real
   kind = kind_phys
   intent = in
-[tskin]
-  standard_name = surface_skin_temperature
-  long_name = surface skin temperature
-  units = K
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
-  intent = in
+
 [hf2d]
   standard_name = instantaneous_surface_upward_sensible_heat_flux
   long_name = surface upward sensible heat flux valid for current call
diff --git a/drivers/ccpp/catchem_gocart_wrapper.F90 b/drivers/ccpp/old/catchem_gocart_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_gocart_wrapper.F90
rename to drivers/ccpp/old/catchem_gocart_wrapper.F90
diff --git a/drivers/ccpp/catchem_gocart_wrapper.meta b/drivers/ccpp/old/catchem_gocart_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_gocart_wrapper.meta
rename to drivers/ccpp/old/catchem_gocart_wrapper.meta
diff --git a/drivers/ccpp/catchem_plume_wrapper.F90 b/drivers/ccpp/old/catchem_plume_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_plume_wrapper.F90
rename to drivers/ccpp/old/catchem_plume_wrapper.F90
diff --git a/drivers/ccpp/catchem_plume_wrapper.meta b/drivers/ccpp/old/catchem_plume_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_plume_wrapper.meta
rename to drivers/ccpp/old/catchem_plume_wrapper.meta
diff --git a/drivers/ccpp/catchem_rad_wrapper.F90 b/drivers/ccpp/old/catchem_rad_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_rad_wrapper.F90
rename to drivers/ccpp/old/catchem_rad_wrapper.F90
diff --git a/drivers/ccpp/catchem_rad_wrapper.meta b/drivers/ccpp/old/catchem_rad_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_rad_wrapper.meta
rename to drivers/ccpp/old/catchem_rad_wrapper.meta
diff --git a/drivers/ccpp/catchem_seas_wrapper.F90 b/drivers/ccpp/old/catchem_seas_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_seas_wrapper.F90
rename to drivers/ccpp/old/catchem_seas_wrapper.F90
diff --git a/drivers/ccpp/catchem_seas_wrapper.meta b/drivers/ccpp/old/catchem_seas_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_seas_wrapper.meta
rename to drivers/ccpp/old/catchem_seas_wrapper.meta
diff --git a/drivers/ccpp/catchem_settling_wrapper.F90 b/drivers/ccpp/old/catchem_settling_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_settling_wrapper.F90
rename to drivers/ccpp/old/catchem_settling_wrapper.F90
diff --git a/drivers/ccpp/catchem_settling_wrapper.meta b/drivers/ccpp/old/catchem_settling_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_settling_wrapper.meta
rename to drivers/ccpp/old/catchem_settling_wrapper.meta
diff --git a/drivers/ccpp/catchem_wetdep_wrapper.F90 b/drivers/ccpp/old/catchem_wetdep_wrapper.F90
similarity index 100%
rename from drivers/ccpp/catchem_wetdep_wrapper.F90
rename to drivers/ccpp/old/catchem_wetdep_wrapper.F90
diff --git a/drivers/ccpp/catchem_wetdep_wrapper.meta b/drivers/ccpp/old/catchem_wetdep_wrapper.meta
similarity index 100%
rename from drivers/ccpp/catchem_wetdep_wrapper.meta
rename to drivers/ccpp/old/catchem_wetdep_wrapper.meta
diff --git a/src/api/catchem.F90 b/src/api/catchem.F90
index 47324865..b4992581 100644
--- a/src/api/catchem.F90
+++ b/src/api/catchem.F90
@@ -11,7 +11,6 @@ module CATChem
    ! CATChem States
    !---------------
    use ChemState_Mod,  only: ChemStateType    !< Chemical State
-   use GridState_Mod,  only: GridStateType    !< Grid State
    use DiagState_Mod,  only: DiagStateType    !< Diagnostic State
    use EmisState_Mod,  only: EmisStateType    !< Emission State
    use MetState_Mod,   only: MetStateType     !< Meteorology State
@@ -50,6 +49,8 @@ module CATChem
    use Error_Mod, only: CC_FAILURE                      !< CATCHem Failure return code
    use Error_Mod, only: CC_SUCCESS                      !< CATChem Successful return code
    use Error_Mod, only: cc_emit_warning => CC_Warning   !< Method for emitting warnings
+   use Error_Mod, only: cc_checkdeallocate => CC_CheckDeallocate !< Method for checking deallocation
+   use Error_Mod, only: cc_checkallocate => CC_CheckAllocate !< Method for checking allocation
    !
    use init_mod, only: cc_init_diag => Init_Diag        !< Method for initializing the diag state
    use init_mod, only: cc_init_met => Init_Met          !< Method for initializing the met state
diff --git a/src/core/chemstate_mod.F90 b/src/core/chemstate_mod.F90
index 00ac1318..e3d857d7 100644
--- a/src/core/chemstate_mod.F90
+++ b/src/core/chemstate_mod.F90
@@ -82,6 +82,54 @@ module ChemState_Mod
 
    end type ChemStateType
 
+   !> \brief Data Type for catchem species
+   !!
+   !! This container defines the catchem species properties
+   !!
+   !!
+   !!!>
+   type, public :: SpeciesType
+
+      ! Names
+      character(len=30) :: long_name  !< long name for species used for netcdf attribute "long_name"
+      character(len=30) :: short_name !< short name for species
+      character(len=50) :: description !< description of species
+
+      ! Logcial switches
+      logical :: is_gas               !< if true, species is a gas and not an aerosol
+      logical :: is_aerosol           !< if true, species is aerosol and not a gas
+      logical :: is_tracer            !< if true, species is a tracer and not an aerosol or gas that undergoes chemistry or photolysis
+      logical :: is_advected          !< if true, species is advected
+      logical :: is_drydep            !< if true, species undergoes dry depotiion
+      logical :: is_photolysis        !< if true, species undergoes photolysis
+      logical :: is_gocart_aero       !< if true, species is a GOCART aerosol species
+      logical :: is_dust              !< if true, species is a dust
+      logical :: is_seasalt           !< if true, species is a seasalt
+
+      ! Numerical properties
+      real(kind=fp) :: mw_g                 !< gaseous molecular weight
+      real(kind=fp) :: density              !< particle density (kg/m3)
+      real(kind=fp) :: radius               !< mean molecular diameter in meters
+      real(kind=fp) :: lower_radius         !< lower radius in meters
+      real(kind=fp) :: upper_radius         !< upper radius in meters
+      real(kind=fp) :: viscosity            !< kinematic viscosity (m2/s)
+
+
+      ! Default background concentration
+      real(kind=fp) :: BackgroundVV        !< Background conc [v/v]
+
+      ! Indices
+      integer :: species_index        !< species index in species array
+      integer :: drydep_index         !< drydep index in drydep array
+      integer :: photolysis_index     !< photolysis index in photolysis array
+      integer :: gocart_aero_index    !< gocart_aero index in gocart_aero array
+
+      ! Concentration
+      real(kind=fp), ALLOCATABLE :: conc(:)             !< species concentration [v/v] or kg/kg
+
+   end type SpeciesType
+
+
 CONTAINS
 
 
@@ -91,24 +139,20 @@ module ChemState_Mod
    !!
    !! \ingroup core_modules
    !!
-   !! \param GridState Grid State
    !! \param ChemState Chem State
    !! \param RC Return code
    !!
    !!!>
-   subroutine Chem_Allocate(GridState, ChemState, RC)
+   subroutine Chem_Allocate(MetState, ChemState, RC)
 
       ! USES
-      USE GridState_Mod,  ONLY : GridStateType
-      USE Species_Mod,    Only : SpeciesType
+      USE MetState_Mod,  ONLY : MetStateType
 
       IMPLICIT NONE
 
       ! INOUT Params
-      type(GridStateType), INTENT(in)    :: GridState ! Grid State object
-      !type(MetStateType), INTENT(in)    :: MetState   ! Met State object
+      type(MetStateType), INTENT(in)    :: MetState   ! Met State object
       type(ChemStateType), INTENT(inout) :: ChemState ! chem State object
-      ! type(SpeciesType),   POINTER       :: Species   ! Species object
       ! OUTPUT Params
       INTEGER,             INTENT(OUT)   :: RC            ! Success or failure
 
@@ -134,7 +178,7 @@ subroutine Chem_Allocate(GridState, ChemState, RC)
       IF ( RC /= CC_SUCCESS ) RETURN
 
       do i=0, ChemState%nSpecies
-         ALLOCATE(ChemState%ChemSpecies(i)%conc(GridState%number_of_levels), STAT=RC)
+         ALLOCATE(ChemState%ChemSpecies(i)%conc(MetState%nLEVS), STAT=RC)
          IF ( RC /= CC_SUCCESS ) THEN
             ErrMsg = 'Could not Allocate ChemState%ChemSpecies(i)%conc'
             CALL CC_Error( ErrMsg, RC, thisLoc )
@@ -144,6 +188,79 @@ subroutine Chem_Allocate(GridState, ChemState, RC)
 
    end subroutine Chem_Allocate
 
+   !> \brief Finalize and deallocate the chem state
+   !!
+   !! \details Deallocates all memory associated with the ChemState object
+   !!
+   !! \ingroup core_modules
+   !!
+   !! \param ChemState Chem State to be deallocated
+   !! \param RC Return code
+   !!
+   !!!>
+   subroutine Chem_Finalize(ChemState, RC)
+      USE CC_Mod,    ONLY : CC_SUCCESS, CC_FAILURE
+      USE CC_Error,  ONLY : CC_CheckDeallocate
+
+      IMPLICIT NONE
+
+      ! Parameters
+      TYPE(ChemStateType), INTENT(INOUT) :: ChemState
+      INTEGER, INTENT(OUT) :: RC
+
+      ! Local variables
+      INTEGER :: i
+
+      ! Initialize
+      RC = CC_SUCCESS
+
+      ! Deallocate concentration arrays for each species
+      do i = 1, ChemState%nSpecies
+         RC = CC_CheckDeallocate(ChemState%ChemSpecies(i)%conc, 'ChemSpecies concentration')
+         if (RC /= CC_SUCCESS) return
+      end do
+
+      ! Deallocate index arrays
+      RC = CC_CheckDeallocate(ChemState%SpeciesIndex, 'SpeciesIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%TracerIndex, 'TracerIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%AeroIndex, 'AeroIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%GasIndex, 'GasIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%DustIndex, 'DustIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%SeaSaltIndex, 'SeaSaltIndex')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(ChemState%DryDepIndex, 'DryDepIndex')
+      if (RC /= CC_SUCCESS) return
+
+      ! Deallocate species names
+      RC = CC_CheckDeallocate(ChemState%SpeciesNames, 'SpeciesNames')
+      if (RC /= CC_SUCCESS) return
+
+      ! Finally deallocate the ChemSpecies array
+      RC = CC_CheckDeallocate(ChemState%ChemSpecies, 'ChemSpecies array')
+      if (RC /= CC_SUCCESS) return
+
+      ! Reset counters to zero
+      ChemState%nSpecies = 0
+      ChemState%nSpeciesGas = 0
+      ChemState%nSpeciesAero = 0
+      ChemState%nSpeciesAeroDryDep = 0
+      ChemState%nSpeciesTracer = 0
+      ChemState%nSpeciesDust = 0
+      ChemState%nSpeciesSeaSalt = 0
+
+   end subroutine Chem_Finalize
+
    !> \brief Find the number of species
    !!
    !! This subroutine finds the number of species
diff --git a/src/core/config_mod.F90 b/src/core/config_mod.F90
index d66b1968..f4da3c64 100644
--- a/src/core/config_mod.F90
+++ b/src/core/config_mod.F90
@@ -38,13 +38,12 @@ MODULE Config_Mod
    !!
    !! \ingroup core_modules
    !!!>
-   SUBROUTINE Read_Input_File( Config , GridState, EmisState, ChemState, RC, ConfigFilename )
+   SUBROUTINE Read_Input_File( Config, EmisState, ChemState, RC, ConfigFilename )
 !
 ! !USES:
 !
       USE Error_Mod
       USE Config_Opt_Mod,  ONLY : ConfigType
-      USE GridState_Mod, ONLY : GridStateType
       use ChemState_Mod, only : ChemStateType
       use EmisState_Mod, only : EmisStateType
 
@@ -54,7 +53,6 @@ SUBROUTINE Read_Input_File( Config , GridState, EmisState, ChemState, RC, Config
 ! !INPUT/OUTPUT PARAMETERS:
 !
       TYPE(ConfigType),    INTENT(INOUT) :: Config    ! Input options
-      TYPE(GridStateType), INTENT(INOUT) :: GridState  ! Grid State object
       TYPE(ChemStateType), INTENT(inout) :: ChemState ! Chemical State
       TYPE(EmisStateType), INTENT(inout) :: EmisState ! Emission State
 !
@@ -120,23 +118,23 @@ SUBROUTINE Read_Input_File( Config , GridState, EmisState, ChemState, RC, Config
          RETURN
       ENDIF
 
-      !========================================================================
-      ! Get grid settings from the YAML Config object
-      !========================================================================
-
-      ! Grid config settings
-      CALL Config_Grid( ConfigInput, GridState, RC )
-      IF ( RC /= CC_SUCCESS ) THEN
-         errMsg = 'Error in "Config_Grid"!'
-         CALL CC_Error( errMsg, RC, thisLoc  )
-         CALL QFYAML_CleanUp( ConfigInput         )
-         CALL QFYAML_CleanUp( ConfigAnchored )
-         RETURN
-      ENDIF
-
       ! !========================================================================
-      ! ! Config processes
+      ! ! Get grid settings from the YAML Config object
       ! !========================================================================
+
+      ! ! Grid config settings
+      ! CALL Config_Grid( ConfigInput, GridState, RC )
+      ! IF ( RC /= CC_SUCCESS ) THEN
+      !    errMsg = 'Error in "Config_Grid"!'
+      !    CALL CC_Error( errMsg, RC, thisLoc  )
+      !    CALL QFYAML_CleanUp( ConfigInput         )
+      !    CALL QFYAML_CleanUp( ConfigAnchored )
+      !    RETURN
+      ! ENDIF
+
+      !========================================================================
+      ! Config processes
+      !========================================================================
       call Config_Process_SeaSalt(ConfigInput, Config, RC)
       IF ( RC /= CC_SUCCESS ) THEN
          errMsg = 'Error in "Config_Process_SeaSalt"!'
@@ -214,15 +212,13 @@ END SUBROUTINE Read_Input_File
    !! \param   RC Return code
    !!
    !!!>
-   SUBROUTINE Config_Chem_State( filename, GridState, ChemState, RC )
+   SUBROUTINE Config_Chem_State( filename, ChemState, RC )
       USE ChemState_Mod, ONLY : ChemStateType, Find_Number_of_Species, Find_Index_of_Species
       use Config_Opt_Mod, ONLY : ConfigType
       USE Error_Mod
-      USE GridState_Mod, ONLY : GridStateType
 
       CHARACTER(LEN=*), INTENT(IN) :: filename
       TYPE(ChemStateType), INTENT(INOUT) :: ChemState
-      TYPE(GridStateType), INTENT(IN) :: GridState
       INTEGER, INTENT(INOUT) :: RC
 
       TYPE(QFYAML_t)     :: ConfigInput, ConfigAnchored
@@ -539,10 +535,10 @@ SUBROUTINE Config_Chem_State( filename, GridState, ChemState, RC )
          ChemState%ChemSpecies(n)%viscosity = v_real
          write(*,*) '|  viscosity: ', ChemState%ChemSpecies(n)%viscosity
 
-         !---------------------------------------
-         ! Allocate initial Species Concentration
-         !---------------------------------------
-         ALLOCATE(ChemState%ChemSpecies(n)%conc(GridState%number_of_levels), STAT=RC)
+         ! !---------------------------------------
+         ! ! Allocate initial Species Concentration
+         ! !---------------------------------------
+         ! ALLOCATE(ChemState%ChemSpecies(n)%conc(GridState%number_of_levels), STAT=RC)
 
       enddo ! n
 
@@ -578,7 +574,6 @@ SUBROUTINE Config_Emis_State( filename, EmisState, RC )
       USE EmisState_Mod, ONLY : EmisStateType
       use Config_Opt_Mod, ONLY : ConfigType
       USE Error_Mod
-      USE GridState_Mod, ONLY : GridStateType
 
       CHARACTER(LEN=*), INTENT(IN) :: filename
       ! TYPE(ChemStateType), INTENT(INOUT) :: ChemState
@@ -853,14 +848,8 @@ SUBROUTINE Config_Simulation( ConfigInput, Config, RC )
 ! !LOCAL VARIABLES:
 !
       ! Scalars
-      ! REAL(fp)                     :: JulianDateStart,  JulianDateEnd
 
       ! Strings
-      ! CHARACTER(LEN=6)             :: timeStr
-      ! CHARACTER(LEN=8)             :: dateStr
-      ! CHARACTER(LEN=12)            :: met
-      ! CHARACTER(LEN=15)            :: verboseMsg
-      ! CHARACTER(LEN=24)            :: sim
       CHARACTER(LEN=255)           :: thisLoc
       CHARACTER(LEN=512)           :: errMsg
       CHARACTER(LEN=QFYAML_NamLen) :: key
@@ -915,109 +904,6 @@ SUBROUTINE Config_Simulation( ConfigInput, Config, RC )
 
    END SUBROUTINE Config_Simulation
 
-   !> \brief Process grid configuration
-   !!
-   !! This function processes the grid configuration and performs the necessary actions based on the configuration.
-   !!
-   !! \param[in] ConfigInput The YAML configuration object
-   !! \param[inout] Config The configuration object
-   !! \param[out] RC The return code
-   !!
-   !! \ingroup core_modules
-   !!!>
-   SUBROUTINE Config_Grid( ConfigInput, GridState, RC )
-!
-! !USES:
-!
-      USE CharPak_Mod,    ONLY : StrSplit
-      USE Error_Mod
-      USE Config_Opt_Mod,  ONLY : ConfigType
-      USE GridState_Mod, ONLY : GridStateType
-!
-! !INPUT/OUTPUT PARAMETERS:
-!
-      TYPE(QFYAML_t),      INTENT(INOUT) :: ConfigInput      ! YAML Config object
-      ! TYPE(ConfigType),     INTENT(INOUT) :: Config   ! Input options
-      TYPE(GridStateType), INTENT(INOUT) :: GridState  ! Grid State
-!
-! !OUTPUT PARAMETERS:
-!
-      INTEGER,        INTENT(OUT)   :: RC          ! Success or failure
-!
-! !LOCAL VARIABLES:
-!
-      ! Scalars
-      ! LOGICAL                      :: v_bool
-      INTEGER                      :: v_int
-
-      ! Strings
-      CHARACTER(LEN=255)           :: thisLoc
-      CHARACTER(LEN=512)           :: errMsg
-      CHARACTER(LEN=QFYAML_StrLen) :: key
-
-      !========================================================================
-      ! Config_Grid begins here!
-      !========================================================================
-
-      ! Initialize
-      RC      = CC_SUCCESS
-      errMsg  = ''
-      thisLoc = ' -> at Config_Grid (in CATChem/src/core/input_mod.F90)'
-
-      !------------------------------------------------------------------------
-      ! Level range
-      !------------------------------------------------------------------------
-      key   = "grid%number_of_levels"
-      v_int = MISSING_INT
-      CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_int, "", RC )
-      IF ( RC /= CC_SUCCESS ) THEN
-         errMsg = 'Error parsing ' // TRIM( key ) // '!'
-         CALL CC_Error( errMsg, RC, thisLoc )
-         RETURN
-      ENDIF
-      GridState%number_of_levels = v_int
-
-      !------------------------------------------------------------------------
-      ! number of soil layers range
-      !------------------------------------------------------------------------
-      key   = "grid%number_of_soil_layers"
-      v_int = MISSING_INT
-      CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_int, "", RC )
-      IF ( RC /= CC_SUCCESS ) THEN
-         errMsg = 'Error parsing ' // TRIM( key ) // '!'
-         CALL CC_Error( errMsg, RC, thisLoc )
-         RETURN
-      ENDIF
-      GridState%number_of_soil_layers = v_int
-
-      !------------------------------------------------------------------------
-      ! number of x and y dimensions (nx and ny)
-      !------------------------------------------------------------------------
-      key   = "grid%nx"
-      v_int = MISSING_INT
-      CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_int, "", RC )
-      IF ( RC /= CC_SUCCESS ) THEN
-         errMsg = 'Error parsing ' // TRIM( key ) // '!'
-         CALL CC_Error( errMsg, RC, thisLoc )
-         RETURN
-      ENDIF
-      GridState%NX = v_int
-
-      key   = "grid%ny"
-      v_int = MISSING_INT
-      CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_int, "", RC )
-      IF ( RC /= CC_SUCCESS ) THEN
-         errMsg = 'Error parsing ' // TRIM( key ) // '!'
-         CALL CC_Error( errMsg, RC, thisLoc )
-         RETURN
-      ENDIF
-      GridState%NY = v_int
-
-      ! Return success
-      RC = CC_SUCCESS
-
-   END SUBROUTINE Config_Grid
-
    !> \brief Process dust configuration
    !!
    !! This function processes the dust configuration and performs the necessary actions based on the configuration.
diff --git a/src/core/diagstate_mod.F90 b/src/core/diagstate_mod.F90
index 1317fd03..f7a0d07b 100644
--- a/src/core/diagstate_mod.F90
+++ b/src/core/diagstate_mod.F90
@@ -15,8 +15,9 @@ module DiagState_Mod
    IMPLICIT NONE
    private
 
-   ! PUBLIC :: Zero_DiagState
    PUBLIC :: Diag_Allocate
+   PUBLIC :: Diag_Finalize
+
 
    !> \brief Data type for storing diagnostic state variables
    !!
@@ -128,4 +129,46 @@ subroutine Diag_Allocate(Config, DiagState, ChemState, RC)
 
    end subroutine Diag_Allocate
 
+
+   !> \brief Finalizes and deallocates the DiagState object
+   !!
+   !! \param State DiagState object to finalize
+   !! \param RC Return code
+   !!
+   !!!>
+   subroutine Diag_Finalize(State, RC)
+      USE CC_Mod,   ONLY : CC_SUCCESS
+      USE CC_Error, ONLY : CC_CheckDeallocate
+
+      !----------------------------------------------------------
+      ! Arguments
+      !----------------------------------------------------------
+      TYPE(DiagStateType), INTENT(INOUT) :: State
+      INTEGER, INTENT(OUT) :: RC
+
+      !----------------------------------------------------------
+      ! Initialize
+      !----------------------------------------------------------
+      RC = CC_SUCCESS
+
+      !----------------------------------------------------------
+      ! Deallocate arrays with error checking
+      !----------------------------------------------------------
+      RC = CC_CheckDeallocate(State%AOD550, 'AOD550')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(State%AOD380, 'AOD380')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(State%TOMSAI, 'TOMSAI')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(State%drydep_frequency, 'drydep_frequency')
+      if (RC /= CC_SUCCESS) return
+
+      RC = CC_CheckDeallocate(State%drydep_vel, 'drydep_vel')
+      if (RC /= CC_SUCCESS) return
+
+   end subroutine Diag_Finalize
+
 end module DiagState_Mod
diff --git a/src/core/emisstate_mod.F90 b/src/core/emisstate_mod.F90
index ec5bcdd7..d14f921f 100644
--- a/src/core/emisstate_mod.F90
+++ b/src/core/emisstate_mod.F90
@@ -135,15 +135,15 @@ module EmisState_Mod
    !! \param RC Return code
    !!
    !!!>
-   subroutine Emis_Allocate(GridState, EmisState, RC)
+   subroutine Emis_Allocate(MetState, EmisState, RC)
 
       ! Uses
-      USE GridState_Mod, ONLY : GridStateType
+      USE MetState_Mod, ONLY : MetStateType
 
       IMPLICIT NONE
 
       !Input
-      TYPE(GridStateType), INTENT(IN) :: GridState
+      TYPE(MetStateType), INTENT(IN) :: MetState
 
       ! Input/Output
       TYPE(EmisStateType), INTENT(INOUT) :: EmisState
@@ -173,7 +173,7 @@ subroutine Emis_Allocate(GridState, EmisState, RC)
             do s = 1, EmisState%Cats(c)%nSpecies
                print*, 'Allocating ', EmisState%Cats(c)%Species(s)%name
 
-               ALLOCATE(EmisState%Cats(c)%Species(s)%Flux(GridState%number_of_levels), STAT=RC)
+               ALLOCATE(EmisState%Cats(c)%Species(s)%Flux(MetState%nLEVS), STAT=RC)
                if (RC /= CC_SUCCESS) then
                   ErrMsg = '  Error allocating "EmisState%Cats%Species%Flux"!'
                   call CC_Error(ErrMsg, RC, ThisLoc)
@@ -448,61 +448,70 @@ end subroutine Apply_Emis_to_Chem
    !! \param RC Return code
    !!
    !!!>
-   subroutine EmisState_CleanUp(EmisState, RC)
+   subroutine EmisState_Finalize(EmisState, RC)
+      USE CC_Mod,   ONLY : CC_SUCCESS
+      USE CC_Error, ONLY : CC_CheckDeallocate
 
       TYPE(EmisStateType), INTENT(INOUT) :: EmisState
       INTEGER, INTENT(OUT) :: RC
 
-      character(len=255) :: ErrMsg
-      character(len=255) :: ThisLoc
-
-      integer :: c ! Loop counter for emission Cats
-      integer :: s ! Loop counter for emitted species
+      integer :: c     ! Loop counter for emission Cats
+      integer :: s     ! Loop counter for emitted species
 
       ! Initialize return code
       RC = CC_SUCCESS
 
-      ! Initialize local variables
-      ErrMsg = ''
-      ThisLoc = ' -> at EmisState_CleanUp (in core/emisstate_mod.F90)'
-
       ! Deallocate total variables
-      if (allocated(EmisState%TotEmisNames)) deallocate(EmisState%TotEmisNames)
+      RC = CC_CheckDeallocate(EmisState%TotEmisNames, 'TotEmisNames')
+      if (RC /= CC_SUCCESS) return
+
+      ! Deallocate total species fluxes
       do c = 1, EmisState%nEmisTotal
-         if (allocated(EmisState%TotSpecies(c)%Flux)) deallocate(EmisState%TotSpecies(c)%Flux)
-         if (allocated(EmisState%TotSpecies(c)%Flux)) deallocate(EmisState%TotSpecies(c)%Flux)
+         RC = CC_CheckDeallocate(EmisState%TotSpecies(c)%Flux, 'TotSpecies Flux')
+         if (RC /= CC_SUCCESS) return
       end do
 
       ! Deallocate emission variables in each category
       cats: do c = 1, EmisState%nCats
          species: do s = 1, EmisState%Cats(c)%nSpecies
-            if (allocated(EmisState%Cats(c)%Species(s)%Flux)) &
-               deallocate(EmisState%Cats(c)%Species(s)%Flux)
-            if (allocated(EmisState%Cats(c)%Species(s)%EmisMapIndex)) &
-               deallocate(EmisState%Cats(c)%Species(s)%EmisMapIndex)
-            if (allocated(EmisState%Cats(c)%Species(s)%Scale)) &
-               deallocate(EmisState%Cats(c)%Species(s)%Scale)
-            if (allocated(EmisState%Cats(c)%Species(s)%EmisMapName)) &
-               deallocate(EmisState%Cats(c)%Species(s)%EmisMapName)
-            if (allocated(EmisState%Cats(c)%Species(s)%PlmRiseHgt)) &
-               deallocate(EmisState%Cats(c)%Species(s)%PlmRiseHgt)
-            if (allocated(EmisState%Cats(c)%Species(s)%PlmSrcFlx)) &
-               deallocate(EmisState%Cats(c)%Species(s)%PlmSrcFlx)
-            if (allocated(EmisState%Cats(c)%Species(s)%FRP)) &
-               deallocate(EmisState%Cats(c)%Species(s)%FRP)
-            if (allocated(EmisState%Cats(c)%Species(s)%STKDM)) &
-               deallocate(EmisState%Cats(c)%Species(s)%STKDM)
-            if (allocated(EmisState%Cats(c)%Species(s)%STKHT)) &
-               deallocate(EmisState%Cats(c)%Species(s)%STKHT)
-            if (allocated(EmisState%Cats(c)%Species(s)%STKTK)) &
-               deallocate(EmisState%Cats(c)%Species(s)%STKTK)
-            if (allocated(EmisState%Cats(c)%Species(s)%STKVE)) &
-               deallocate(EmisState%Cats(c)%Species(s)%STKVE)
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%Flux, 'Species Flux')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%EmisMapIndex, 'EmisMapIndex')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%Scale, 'Scale')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%EmisMapName, 'EmisMapName')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%PlmRiseHgt, 'PlmRiseHgt')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%PlmSrcFlx, 'PlmSrcFlx')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%FRP, 'FRP')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%STKDM, 'STKDM')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%STKHT, 'STKHT')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%STKTK, 'STKTK')
+            if (RC /= CC_SUCCESS) return
+
+            RC = CC_CheckDeallocate(EmisState%Cats(c)%Species(s)%STKVE, 'STKVE')
+            if (RC /= CC_SUCCESS) return
          end do species
       end do cats
-      if (allocated(EmisState%Cats)) deallocate(EmisState%Cats)
-      if (allocated(EmisState%Cats)) deallocate(EmisState%Cats)
 
-   end subroutine EmisState_CleanUp
+      ! Final deallocation of Cats array
+      RC = CC_CheckDeallocate(EmisState%Cats, 'Cats array')
+      if (RC /= CC_SUCCESS) return
 
+   end subroutine EmisState_Finalize
 END MODULE EmisState_Mod
diff --git a/src/core/error_mod.F90 b/src/core/error_mod.F90
index e6f32f1f..817a1150 100644
--- a/src/core/error_mod.F90
+++ b/src/core/error_mod.F90
@@ -19,6 +19,7 @@ MODULE Error_Mod
    PUBLIC :: CC_Error
    PUBLIC :: CC_Warning
    PUBLIC :: CC_CheckVar
+   PUBLIC :: CC_CheckDeallocate
    !
    ! !DEFINED PARAMETERS:
    !
@@ -212,5 +213,102 @@ SUBROUTINE CC_CheckVar( Variable, Operation, RC )
       ENDIF
 
    END SUBROUTINE CC_CheckVar
+
+   !<
+   !! \brief Check and perform array allocation
+   !!
+   !! \param array    Array to allocate
+   !! \param size     Size of array to allocate
+   !! \param varName  Variable name for error messages
+   !!
+   !! \return RC     Return code (CC_SUCCESS or error)
+   !!!>
+   FUNCTION CC_CheckAllocate(array, size, varName) RESULT(RC)
+      !
+      ! !INPUT PARAMETERS:
+      !
+      CLASS(*), ALLOCATABLE, INTENT(INOUT) :: array(:)  !< Array to allocate
+      INTEGER,              INTENT(IN)    :: size      !< Size to allocate
+      CHARACTER(LEN=*),     INTENT(IN)    :: varName   !< Variable name for messages
+      !
+      ! !RETURN VALUE:
+      !
+      INTEGER :: RC
+      !
+      ! !LOCAL VARIABLES:
+      !
+      INTEGER           :: stat
+      CHARACTER(LEN=255) :: ErrMsg, ThisLoc
+
+      !=========================================================================
+      ! Initialize
+      !=========================================================================
+      RC = CC_SUCCESS
+      stat = 0
+      ThisLoc = ' -> at CC_CheckAllocate (in Headers/error_mod.F90)'
+
+      !=========================================================================
+      ! Perform allocation if array is not already allocated
+      !=========================================================================
+      IF (.NOT. ALLOCATED(array)) THEN
+         ALLOCATE(array(size), STAT=stat)
+
+         IF (stat /= 0) THEN
+            ErrMsg = 'Failed to allocate ' // TRIM(varName)
+            CALL CC_Error(ErrMsg, RC, ThisLoc)
+            RETURN
+         ENDIF
+      ENDIF
+
+   END FUNCTION CC_CheckAllocate
+
+   !>
+   !! \brief CC_CheckDeallocate
+   !!
+   !! This function safely deallocates arrays and returns appropriate error codes
+   !!
+   !! \param array Array to deallocate
+   !! \param varName Name of the variable being deallocated (for error messages)
+   !! \return RC Return code (CC_SUCCESS or CC_FAILURE)
+   !!
+   !! \ingroup core_modules
+   !!!>
+   FUNCTION CC_CheckDeallocate(array, varName) RESULT(RC)
+      !
+      ! !INPUT PARAMETERS:
+      !
+      CLASS(*), ALLOCATABLE, INTENT(INOUT) :: array    !< Array to deallocate
+      CHARACTER(LEN=*),      INTENT(IN)    :: varName  !< Variable name for messages
+      !
+      ! !RETURN VALUE:
+      !
+      INTEGER :: RC
+      !
+      ! !LOCAL VARIABLES:
+      !
+      INTEGER          :: stat
+      CHARACTER(LEN=255) :: ErrMsg, ThisLoc
+
+      !=========================================================================
+      ! Initialize
+      !=========================================================================
+      RC = CC_SUCCESS
+      stat = 0
+      ThisLoc = ' -> at CC_CheckDeallocate (in Headers/error_mod.F90)'
+
+      !=========================================================================
+      ! Perform deallocation if array is allocated
+      !=========================================================================
+      IF (ALLOCATED(array)) THEN
+         DEALLOCATE(array, STAT=stat)
+
+         IF (stat /= 0) THEN
+            ErrMsg = 'Failed to deallocate ' // TRIM(varName)
+            CALL CC_Error(ErrMsg, RC, ThisLoc)
+            RETURN
+         ENDIF
+      ENDIF
+
+   END FUNCTION CC_CheckDeallocate
    !EOC
 END MODULE Error_Mod
diff --git a/src/core/gridstate_mod.F90 b/src/core/gridstate_mod.F90
deleted file mode 100644
index c2b9e91f..00000000
--- a/src/core/gridstate_mod.F90
+++ /dev/null
@@ -1,70 +0,0 @@
-!> \file gridstate_mod.F90
-!!
-!! \brief Module for grid state variables
-!!
-!! This module contains subroutines and functions related to the grid state.
-!!
-!! \ingroup core_modules
-!!!>
-module GridState_Mod
-
-   USE Error_Mod
-   USE precision_mod
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC :: Grid_Init_State
-   type, public :: GridStateType
-      CHARACTER(LEN=4) :: State = 'Grid'  !< Name of this state
-
-      ! Integers
-      integer :: nx = 1
-      integer :: ny = 1
-      integer :: number_of_levels  !< The number of vertical levels
-      integer :: number_of_soil_layers  !< The number of soil layers
-
-      ! Reals
-      real(fp) :: area  !< Grid cell horizontal area [m^2]
-
-   end type GridStateType
-
-contains
-
-   !> \brief Initialize a GridState object
-   !!
-   !! This subroutine initializes a GridState object.
-   !!
-   !! \param Config The input config object.
-   !! \param GridState The GridState object to be initialized.
-   !! \param RC The return code.
-   !!
-   !! \ingroup core_modules
-   !!!>
-   subroutine Grid_Init_State(GridState, RC)
-      use Error_Mod, only : CC_SUCCESS
-      use Config_Opt_Mod, Only : ConfigType
-      implicit none
-
-      ! type(ConfigType),    intent(in)    :: Config     !< Input Options object
-      type(GridStateType), intent(inout) :: GridState  !< Grid State object
-      INTEGER,             INTENT(OUT)   :: RC         !< Success or failure
-
-      ! Local variables
-      CHARACTER(LEN=512) :: errMsg
-      CHARACTER(LEN=255) :: thisLoc
-
-      ! Set error handling defaults
-      RC = CC_SUCCESS
-      errMsg = ''
-      thisLoc = 'Grid_Init_State() -> at initializing GridState'
-
-      ! initialize GridState
-      GridState%nx=1
-      GridState%ny=1
-      GridState%number_of_levels=1  ! FIXME: use Config?
-      GridState%area = 1._fp
-
-   end subroutine Grid_Init_State
-
-end module GridState_Mod
diff --git a/src/core/metstate_mod.F90 b/src/core/metstate_mod.F90
index 164ccdb1..cbe4ecb4 100644
--- a/src/core/metstate_mod.F90
+++ b/src/core/metstate_mod.F90
@@ -10,11 +10,8 @@ MODULE MetState_Mod
    !
    ! USES:
    !
-   USE Cmn_Size_Mod, ONLY : NSURFTYPE
-   ! USE Dictionary_M, ONLY : dictionary_t
    USE Error_Mod
    USE Precision_Mod
-   ! USE Registry_Mod
 
 
    IMPLICIT NONE
@@ -38,25 +35,36 @@ MODULE MetState_Mod
 
       CHARACTER(LEN=3)             :: State     = 'MET'    ! Name of this state
 
-      ! NLEVS
-      !------
-      INTEGER               :: NLEVS             !< Number of vertical levels
+      ! Integer Fields for MetState Array Dimensions
+      !---------------------------------------------
+      INTEGER               :: nLEVS             !< Number of vertical levels
+      INTEGER               :: nSOIL             !< # number of soil layers
+      INTEGER               :: nLNDTYPE          !< # of landtypes in box (I,J)
+
+      ! Location Specific Fields
+      !-------------------------
+      real(fp) :: LAT               !< Latitude [degrees]
+      real(fp) :: LON               !< Longitude [degrees]
 
       ! TIMESTEP
       !---------
-      REAL(fp)              :: TSTEP             !< Time step [s]
+      REAL(fp) :: TSTEP             !< Time step [s]
+      REAL(fp) :: JDAY              !< Julian Day of year
+      REAL(fp) :: JDAY_FRAC         !< Fractional Julian Day of year
+      real(fp) :: HR                !< Hour of Day in UTC
 
       ! Logicals
       !---------
-      LOGICAL           :: IsLand            !< Is this a land grid box?
-      LOGICAL           :: IsWater           !< Is this a water grid box?
-      LOGICAL           :: IsIce             !< Is this a ice grid box?
-      LOGICAL           :: IsSnow            !< Is this a snow grid box?
+      LOGICAL           :: IsLand                !< Is this a land grid box?
+      LOGICAL           :: IsWater               !< Is this a water grid box?
+      LOGICAL           :: IsIce                 !< Is this a ice grid box?
+      LOGICAL           :: IsSnow                !< Is this a snow grid box?
+      LOGICAL           :: IsLocalNoon           ! Is it local noon (between 11 and 13 local solar time?
       LOGICAL,  ALLOCATABLE :: InStratMeso(:)    !< Are we in the stratosphere or mesosphere?
       LOGICAL,  ALLOCATABLE :: InStratosphere(:) !< Are we in the stratosphere?
       LOGICAL,  ALLOCATABLE :: InTroposphere(:)  !< Are we in the troposphere?
       LOGICAL,  ALLOCATABLE :: InPbl(:)          !< Are we in the PBL?
-      LOGICAL,  ALLOCATABLE :: IsLocalNoon       ! Is it local noon (between 11 and 13 local solar time?
+
 
       ! Land Specific Fields
       !---------------------
@@ -70,8 +78,9 @@ MODULE MetState_Mod
       REAL(fp)              :: FRLAND          !< Fraction of land [1]
       REAL(fp)              :: FRLANDIC        !< Fraction of land ice [1]
       REAL(fp)              :: FROCEAN         !< Fraction of ocean [1]
-      REAL(fp)              :: FRSEAICE        !< Sfc sea ice fraction
-      REAL(fp)              :: FRSNO           !< Sfc snow fraction
+      REAL(fp)              :: FRSEAICE        !< Sfc sea ice fraction [1]
+      REAL(fp)              :: FRSNO           !< Sfc snow fraction [1]
+      REAL(fp)              :: FRURBAN         !< Fraction of urban [1]
       REAL(fp)              :: LAI             !< Leaf area index [m2/m2] (online) Dominant
       REAL(fp)              :: GVF             !< Green Vegetative Fraction
       REAL(fp)              :: RDRAG           !< Drag Partition [1]
@@ -90,12 +99,11 @@ MODULE MetState_Mod
       REAL(fp)              :: SNOMAS          !< Snow mass [kg/m2]
       REAL(fp)              :: SSM             !< Sediment Supply Map [1]
       REAL(fp)              :: USTAR_THRESHOLD !< Threshold friction velocity [m/s]
-      INTEGER,  ALLOCATABLE :: nLNDTYPE        !< # of landtypes in box (I,J)
       REAL(fp)              :: GWETTOP         !< Top soil moisture [1]
       REAL(fp)              :: GWETROOT        !< Root Zone soil moisture [1]
       REAL(fp)              :: WILT            !< Wilt point [1]
-      INTEGER,  ALLOCATABLE :: nSOIL           !< # number of soil layers
       REAL(fp), ALLOCATABLE :: SOILM(:)        !< Volumetric Soil moisture [m3/m3]
+      REAL(fp), ALLOCATABLE :: SOILT(:)        !< Volumetric Soil T [K]
       REAL(fp), ALLOCATABLE :: FRLANDUSE(:)    !< Fractional Land Use
       REAL(fp), ALLOCATABLE :: FRLAI(:)        !< LAI in each Fractional Land use type [m2/m2]
 
@@ -224,24 +232,13 @@ SUBROUTINE Zero_MetState( MetState, RC )
 
    END SUBROUTINE Zero_MetState
 
-   !>
-   !! \brief Allocate the MetState object
-   !!
-   !! \ingroup core_modules
-   !!
-   !! \param GridState   CATCHem grid state
-   !! \param MetState    CATCHem met state
-   !! \param RC          Error return code
-   !!!>
-   SUBROUTINE Met_Allocate( GridState, MetState, RC)
-      ! USES
-      USE GridState_Mod, Only : GridStateType
+   SUBROUTINE Met_Allocate(MetState, RC)
+      USE error_mod,     ONLY : CC_CheckAllocate, CC_SUCCESS, CC_Error
 
       IMPLICIT NONE
 
       ! Arguments
-      TYPE(GridStateType), INTENT(IN)  :: GridState !< Grid state
-      TYPE(MetStateType), INTENT(INOUT) :: MetState !< Meteorological state
+      TYPE(MetStateType),  INTENT(INOUT) :: MetState !< Meteorological state
       INTEGER,            INTENT(OUT)   :: RC       !< Return code
 
       ! Local variables
@@ -252,16 +249,6 @@ SUBROUTINE Met_Allocate( GridState, MetState, RC)
       ErrMsg = ''
       thisLoc = ' -> at Met_Allocate (in core/metstate_mod.F90)'
 
-      ! Nullify all fields for safety's sake before allocating them
-      ! This can prevent compilation errors caused by uninitialized values
-
-
-      !--------------------------------------------------
-      ! Initialize fields
-      !--------------------------------------------------
-      MetState%nSOIL = GridState%number_of_soil_layers
-      print*, 'MetState%nSOIL = ', MetState%nSOIL
-
       ! Visible Surface Albedo
       !-----------------------
       MetState%ALBD_VIS = ZERO
@@ -303,400 +290,275 @@ SUBROUTINE Met_Allocate( GridState, MetState, RC)
       MetSTate%SANDFRAC = ZERO
       MetState%SST = ZERO
 
-      ! Allocate Column Fields
+      ! Allocate Column Fields using CC_CheckAllocate
       !-----------------------
-      !  Logicals
-      if (.not. allocated(MetState%InStratosphere)) then
-         allocate(MetState%InStratosphere(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%InPbl)) then
-         allocate(MetState%InPbl(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InPbl'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%InStratMeso)) then
-         allocate(MetState%InStratMeso(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratMeso'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%InTroposphere)) then
-         allocate(MetState%InTroposphere(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InTroposphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
+      ! Logicals
+      RC = CC_CheckAllocate(MetState%InStratosphere, MetState%nLEVS, 'MetState%InStratosphere')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%InPbl, MetState%nLEVS, 'MetState%InPbl')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%InStratMeso, MetState%nLEVS, 'MetState%InStratMeso')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%InTroposphere, MetState%nLEVS, 'MetState%InTroposphere')
+      IF (RC /= CC_SUCCESS) RETURN
 
       ! Flux Related
-      if (.not. allocated(MetState%F_OF_PBL)) then
-         allocate(MetState%F_OF_PBL(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%F_OF_PBL'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%F_UNDER_PBLTOP)) then
-         allocate(MetState%F_UNDER_PBLTOP(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%F_UNDER_PBLTOP'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
+      RC = CC_CheckAllocate(MetState%F_OF_PBL, MetState%nLEVS, 'MetState%F_OF_PBL')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%F_UNDER_PBLTOP, MetState%nLEVS, 'MetState%F_UNDER_PBLTOP')
+      IF (RC /= CC_SUCCESS) RETURN
 
       ! Cloud / Precipitation
-      if (.not. allocated(MetState%CLDF)) then
-         allocate(MetState%CLDF(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%CLDF'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%CMFMC)) then
-         allocate(MetState%CMFMC(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%CMFMC'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DQRCU)) then
-         allocate(MetState%DQRCU(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%DQRCU'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DQRLSAN)) then
-         allocate(MetState%DQRLSAN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%DQRLSAN'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DTRAIN)) then
-         allocate(MetState%DTRAIN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%DTRAIN'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%QI)) then
-         allocate(MetState%QI(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%QI'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%QL)) then
-         allocate(MetState%QL(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%QL'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PFICU)) then
-         allocate(MetState%PFICU(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%PFICU'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PFILSAN)) then
-         allocate(MetState%PFILSAN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%PFILSAN'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PFLCU)) then
-         allocate(MetState%PFLCU(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%PFLCU'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PFLLSAN)) then
-         allocate(MetState%PFLLSAN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%PFLLSAN'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%TAUCLI)) then
-         allocate(MetState%TAUCLI(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%TAUCLI'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%TAUCLW)) then
-         allocate(MetState%TAUCLW(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%TAUCLW'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      ! State Variables
-      if (.not. allocated(MetState%Z)) then
-         allocate(MetState%Z(GridState%number_of_levels + 1), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%Z'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%ZMID)) then
-         allocate(MetState%ZMID(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%ZMID'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%BXHEIGHT)) then
-         allocate(MetState%BXHEIGHT(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%BXHEIGHT'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%QV)) then
-         allocate(MetState%QV(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%T)) then
-         allocate(MetState%T(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%T'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%THETA)) then
-         allocate(MetState%THETA(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%THETA'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%TV)) then
-         allocate(MetState%TV(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%TV'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%U)) then
-         allocate(MetState%U(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%V)) then
-         allocate(MetState%V(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%OMEGA)) then
-         allocate(MetState%OMEGA(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%OMEGA'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%RH)) then
-         allocate(MetState%RH(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%SPHU)) then
-         allocate(MetState%SPHU(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%AIRDEN)) then
-         allocate(MetState%AIRDEN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%AIRNUMDEN)) then
-         allocate(MetState%AIRNUMDEN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%MAIRDEN)) then
-         allocate(MetState%MAIRDEN(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%AVGW)) then
-         allocate(MetState%AVGW(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DELP)) then
-         allocate(MetState%DELP(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DELP_DRY)) then
-         allocate(MetState%DELP_DRY(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%DAIRMASS)) then
-         allocate(MetState%DAIRMASS(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%AIRVOL)) then
-         allocate(MetState%AIRVOL(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PMID)) then
-         allocate(MetState%PMID(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PMID_DRY)) then
-         allocate(MetState%PMID_DRY(GridState%number_of_levels), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%PEDGE_DRY)) then
-         allocate(MetState%PEDGE_DRY(GridState%number_of_levels+1), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-      if (.not. allocated(MetState%SOILM)) then
-         allocate(MetState%SOILM(MetState%nSOIL), stat=RC)
-         if (RC /= CC_SUCCESS) then
-            errMsg = 'Error allocating MetState%InStratosphere'
-            call CC_Error(errMsg, RC, thisLoc)
-            return
-         endif
-      end if
-
-   end subroutine Met_Allocate
+      RC = CC_CheckAllocate(MetState%CLDF, MetState%nLEVS, 'MetState%CLDF')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%CMFMC, MetState%nLEVS, 'MetState%CMFMC')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DQRCU, MetState%nLEVS, 'MetState%DQRCU')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DQRLSAN, MetState%nLEVS, 'MetState%DQRLSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DTRAIN, MetState%nLEVS, 'MetState%DTRAIN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%QI, MetState%nLEVS, 'MetState%QI')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%QL, MetState%nLEVS, 'MetState%QL')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PFICU, MetState%nLEVS, 'MetState%PFICU')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PFILSAN, MetState%nLEVS, 'MetState%PFILSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PFLCU, MetState%nLEVS, 'MetState%PFLCU')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PFLLSAN, MetState%nLEVS, 'MetState%PFLLSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%TAUCLI, MetState%nLEVS, 'MetState%TAUCLI')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%TAUCLW, MetState%nLEVS, 'MetState%TAUCLW')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      ! State Related
+      RC = CC_CheckAllocate(MetState%Z, MetState%nLEVS + 1, 'MetState%Z')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%ZMID, MetState%nLEVS, 'MetState%ZMID')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%BXHEIGHT, MetState%nLEVS, 'MetState%BXHEIGHT')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%QV, MetState%nLEVS, 'MetState%QV')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%T, MetState%nLEVS, 'MetState%T')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%THETA, MetState%nLEVS, 'MetState%THETA')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%TV, MetState%nLEVS, 'MetState%TV')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%U, MetState%nLEVS, 'MetState%U')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%V, MetState%nLEVS, 'MetState%V')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%OMEGA, MetState%nLEVS, 'MetState%OMEGA')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%RH, MetState%nLEVS, 'MetState%RH')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%SPHU, MetState%nLEVS, 'MetState%SPHU')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%AIRDEN, MetState%nLEVS, 'MetState%AIRDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%AIRNUMDEN, MetState%nLEVS, 'MetState%AIRNUMDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%MAIRDEN, MetState%nLEVS, 'MetState%MAIRDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%AVGW, MetState%nLEVS, 'MetState%AVGW')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DELP, MetState%nLEVS, 'MetState%DELP')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DELP_DRY, MetState%nLEVS, 'MetState%DELP_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%DAIRMASS, MetState%nLEVS, 'MetState%DAIRMASS')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%AIRVOL, MetState%nLEVS, 'MetState%AIRVOL')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PMID, MetState%nLEVS, 'MetState%PMID')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PMID_DRY, MetState%nLEVS, 'MetState%PMID_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%PEDGE_DRY, MetState%nLEVS + 1, 'MetState%PEDGE_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      ! Surface and Soil Properties
+      RC = CC_CheckAllocate(MetState%SOILM, MetState%nSOIL, 'MetState%SOILM')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%SOILT, MetState%nSOIL, 'MetState%SOILT')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%FRLANDUSE, MetState%nLNDTYPE, 'MetState%FRLANDUSE')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%FRLAI, MetState%nLNDTYPE, 'MetState%FRLAI')
+      IF (RC /= CC_SUCCESS) RETURN
+
+      RC = CC_CheckAllocate(MetState%FRZ0, MetState%nLNDTYPE, 'MetState%FRZ0')
+      IF (RC /= CC_SUCCESS) RETURN
+
+   END SUBROUTINE Met_Allocate
+
+   !>
+   !! \brief Deallocate the MetState object
+   !!
+   !! \ingroup core_modules
+   !!
+   !! \param MetState    CATCHem met state
+   !! \param RC          Error return code
+   !!!>
+   SUBROUTINE Met_Finalize( MetState, RC )
+      ! Arguments
+      TYPE(MetStateType), INTENT(INOUT) :: MetState !< Meteorological state
+      INTEGER,            INTENT(OUT)   :: RC       !< Return code
+
+      ! Local variables
+      CHARACTER(LEN=255) :: ErrMsg, thisLoc
+
+      ! Initialize
+      RC = CC_SUCCESS
+      ErrMsg = ''
+      thisLoc = ' -> at Met_Finalize (in core/metstate_mod.F90)'
+
+      ! Deallocate all allocated arrays
+      RC = CC_CheckDeallocate(MetState%InStratosphere, 'MetState%InStratosphere')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%InPbl,          'MetState%InPbl')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%InStratMeso,    'MetState%InStratMeso')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%InTroposphere,  'MetState%InTroposphere')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%F_OF_PBL,       'MetState%F_OF_PBL')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%F_UNDER_PBLTOP, 'MetState%F_UNDER_PBLTOP')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%CLDF,           'MetState%CLDF')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%CMFMC,          'MetState%CMFMC')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DQRCU,          'MetState%DQRCU')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DQRLSAN,        'MetState%DQRLSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DTRAIN,         'MetState%DTRAIN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%QI,             'MetState%QI')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%QL,             'MetState%QL')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PFICU,          'MetState%PFICU')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PFILSAN,        'MetState%PFILSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PFLCU,          'MetState%PFLCU')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PFLLSAN,        'MetState%PFLLSAN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%TAUCLI,         'MetState%TAUCLI')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%TAUCLW,         'MetState%TAUCLW')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%Z,              'MetState%Z')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%ZMID,           'MetState%ZMID')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%BXHEIGHT,       'MetState%BXHEIGHT')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%QV,             'MetState%QV')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%T,              'MetState%T')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%THETA,          'MetState%THETA')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%TV,             'MetState%TV')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%U,              'MetState%U')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%V,              'MetState%V')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%OMEGA,          'MetState%OMEGA')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%RH,             'MetState%RH')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%SPHU,           'MetState%SPHU')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%AIRDEN,         'MetState%AIRDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%AIRNUMDEN,      'MetState%AIRNUMDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%MAIRDEN,        'MetState%MAIRDEN')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%AVGW,           'MetState%AVGW')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DELP,           'MetState%DELP')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DELP_DRY,       'MetState%DELP_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%DAIRMASS,       'MetState%DAIRMASS')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%AIRVOL,         'MetState%AIRVOL')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PMID,           'MetState%PMID')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PMID_DRY,       'MetState%PMID_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%PEDGE_DRY,      'MetState%PEDGE_DRY')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%SOILM,          'MetState%SOILM')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%SOILT,          'MetState%SOILT')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%FRLANDUSE,      'MetState%FRLANDUSE')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%FRLAI,          'MetState%FRLAI')
+      IF (RC /= CC_SUCCESS) RETURN
+      RC = CC_CheckDeallocate(MetState%FRZ0,           'MetState%FRZ0')
+      IF (RC /= CC_SUCCESS) RETURN
+
+   END SUBROUTINE Met_Finalize
+
+
 
 END MODULE MetState_Mod
diff --git a/src/core/species_mod.F90 b/src/core/species_mod.F90
deleted file mode 100644
index 338d18bc..00000000
--- a/src/core/species_mod.F90
+++ /dev/null
@@ -1,101 +0,0 @@
-!>
-!! \file species_mod.F90
-!! \brief This file contains the module for catchem species
-!!
-!! \ingroup core_modules
-!!
-!! This file contains the module for catchem species
-!!
-!!!>
-
-module species_mod
-
-   use precision_mod
-   implicit none
-
-   !> \brief Module for catchem species
-   !!
-   !! This module contains subroutines and functions related to the catchem species
-   !!
-   !! \ingroup core_modules
-   !!
-   !! \param Config The input config object.
-   !! \param Species The Species object to be initialized.
-   !! \param RC The return code.
-   !!
-   !!!>
-   type, public :: SpeciesType
-
-      ! Names
-      character(len=30) :: long_name  !< long name for species used for netcdf attribute "long_name"
-      character(len=30) :: short_name !< short name for species
-      character(len=50) :: description !< description of species
-
-      ! Logcial switches
-      logical :: is_gas               !< if true, species is a gas and not an aerosol
-      logical :: is_aerosol           !< if true, species is aerosol and not a gas
-      logical :: is_tracer            !< if true, species is a tracer and not an aerosol or gas that undergoes chemistry or photolysis
-      logical :: is_advected          !< if true, species is advected
-      logical :: is_drydep            !< if true, species undergoes dry depotiion
-      logical :: is_photolysis        !< if true, species undergoes photolysis
-      logical :: is_gocart_aero       !< if true, species is a GOCART aerosol species
-      logical :: is_dust              !< if true, species is a dust
-      logical :: is_seasalt           !< if true, species is a seasalt
-
-      ! Numerical properties
-      real(kind=fp) :: mw_g                 !< gaseous molecular weight
-      real(kind=fp) :: density              !< particle density (kg/m3)
-      real(kind=fp) :: radius               !< mean molecular diameter in meters
-      real(kind=fp) :: lower_radius         !< lower radius in meters
-      real(kind=fp) :: upper_radius         !< upper radius in meters
-      real(kind=fp) :: viscosity            !< kinematic viscosity (m2/s)
-
-
-      ! Default background concentration
-      real(kind=fp) :: BackgroundVV        !< Background conc [v/v]
-
-      ! Indices
-      integer :: species_index        !< species index in species array
-      integer :: drydep_index         !< drydep index in drydep array
-      integer :: photolysis_index     !< photolysis index in photolysis array
-      integer :: gocart_aero_index    !< gocart_aero index in gocart_aero array
-
-      ! Concentration
-      real(kind=fp), ALLOCATABLE :: conc(:)             !< species concentration [v/v] or kg/kg
-
-   end type SpeciesType
-
-   !
-   ! !DEFINED PARAMETERS:
-   !
-   !=========================================================================
-   ! Missing species concentration value if not in restart file and special
-   ! background value not defined
-   !=========================================================================
-   REAL(fp), PARAMETER, PUBLIC :: MISSING_VV  = 1.0e-20_fp ! Missing spc conc
-
-contains
-
-   subroutine init(Species_State, species_name, atomic_num)
-      type(SpeciesType), intent(inout) :: Species_State
-      character(len=*), intent(in) :: species_name
-      integer, intent(in) :: atomic_num
-
-      Species_State%short_name = species_name
-      Species_State%mw_g = atomic_num
-   end subroutine init
-
-   ! function get_name(this) result(species_name)
-   !    character(len=30) :: species_name
-
-   !    species_name = this%short_name
-   ! end function get_name
-
-   ! function get_atomic_number(this) result(atomic_num)
-   !    class(Species), intent(in) :: this
-   !    integer :: atomic_num
-
-   !    atomic_num = this%atomic_number
-   ! end function get_atomic_number
-
-end module species_mod
diff --git a/src/core/state_mod.F90 b/src/core/state_mod.F90
index d591ee98..52879bb4 100644
--- a/src/core/state_mod.F90
+++ b/src/core/state_mod.F90
@@ -9,7 +9,6 @@
 module state_mod
    use precision_mod
    use Config_Opt_Mod, only : ConfigType
-   use GridState_Mod,  only : GridStateType
    use MetState_Mod,   only : MetStateType
    use ChemState_Mod,  only : ChemStateType
    use EmisState_Mod,  only : EmisStateType
@@ -18,7 +17,6 @@ module state_mod
    IMPLICIT NONE
 
    ! PUBLIC
-   type(GridStateType), PUBLIC :: GridState
    type(MetStateType),  PUBLIC :: MetState
    type(ChemStateType), PUBLIC :: ChemState
    type(ConfigType),    PUBLIC :: Config
diff --git a/tests/test_dust.f90 b/tests/test_dust.f90
index 6c07ef5f..efbf8662 100644
--- a/tests/test_dust.f90
+++ b/tests/test_dust.f90
@@ -35,7 +35,7 @@ program test_dust
    !----------------------------
 
    ! Read input file and initialize grid
-   call cc_read_config(Config, GridState, EmisState, ChemState, rc, configFile)
+   call cc_read_config(Config, EmisState, ChemState, rc, configFile)
    if (rc /= CC_success) then
       errMsg = 'Error reading configuration file: ' // TRIM( configFile )
       call cc_emit_error(errMsg, rc, thisLoc)
diff --git a/tests/test_main.F90 b/tests/test_main.F90
index 25ed22de..a1bd6618 100644
--- a/tests/test_main.F90
+++ b/tests/test_main.F90
@@ -36,7 +36,7 @@ program test_main
    write(*,*) ''
 
    ! Read input file and initialize grid
-   call cc_read_config(Config, GridState, EmisState, ChemState, rc, configFile)
+   call cc_read_config(Config, EmisState, ChemState, rc, configFile)
    if (rc /= CC_SUCCESS) then
       errMsg = 'Error reading configuration file: ' // TRIM( configFile )
       call cc_emit_error(errMsg, rc, thisLoc)
@@ -99,14 +99,8 @@ program test_main
    call assert_close(ChemState%ChemSpecies(1)%radius, 0.8_fp, msg="dust1 radius")
    call assert_close(ChemState%ChemSpecies(1)%density, 2500.0_fp, msg="dust1 density")
 
-   ! write grid info
-   write(*,*) 'Grid info:'
-   write(*,*) 'Number of grid nx = ', GridState%NX
-   write(*,*) 'Number of grid ny = ', GridState%NY
-   write(*,*) 'Number of grid levels = ', GridState%number_of_levels
-
    ! initialize met
-   call cc_init_met(GridState, MetState, rc)
+   call cc_init_met(MetState, rc)
    if (rc /= CC_SUCCESS) then
       errMsg = 'Error initializing meteorology'
       call cc_emit_error(errMsg, rc, thisLoc)
diff --git a/tests/test_plumerise.F90 b/tests/test_plumerise.F90
index 163a370d..72c8b22d 100644
--- a/tests/test_plumerise.F90
+++ b/tests/test_plumerise.F90
@@ -35,7 +35,7 @@ program test_plumerise
    write(*,*) '  PLUMERISE TEST'
 
    ! Read input file and initialize grid
-   call cc_read_config(Config, GridState, EmisState, ChemState, rc, configFile)
+   call cc_read_config(Config, EmisState, ChemState, rc, configFile)
    if (rc /= CC_success) then
       errMsg = 'Error reading configuration file: ' // TRIM( configFile )
       call cc_emit_error(errMsg, rc, thisLoc)
@@ -44,8 +44,9 @@ program test_plumerise
 
    ! Allocate MetState
    MetState%nSOIL = 4
+   MetState%nLEVS = 127
    print*, 'Allocating MetState'
-   call cc_allocate_metstate(GridState, MetState, rc)
+   call cc_allocate_metstate(MetState, rc)
    if (rc /= CC_success) then
       errMsg = 'Error in "cc_allocate_metstate"'
       call cc_emit_error(errMsg, rc, thisLoc)
@@ -62,7 +63,7 @@ program test_plumerise
    endif
 
    ! Allocate EmisState
-   call cc_allocate_emisstate(GridState, EmisState, rc)
+   call cc_allocate_emisstate(MetState, EmisState, rc)
    if (rc /= CC_success) then
       errMsg = 'Error in "cc_allocate_emisstate"'
       call cc_emit_error(errMsg, rc, thisLoc)
diff --git a/tests/test_seasalt.F90 b/tests/test_seasalt.F90
index ec16da80..82e12707 100644
--- a/tests/test_seasalt.F90
+++ b/tests/test_seasalt.F90
@@ -35,7 +35,7 @@ program test_dust
    !----------------------------
 
    ! Read input file and initialize grid
-   call cc_read_config(Config, GridState, EmisState, ChemState, rc, configFile)
+   call cc_read_config(Config, EmisState, ChemState, rc, configFile)
    if (rc /= CC_success) then
       errMsg = 'Error reading configuration file: ' // TRIM( configFile )
       call cc_emit_error(errMsg, rc, thisLoc)