diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 9663a3cbf3..d25920e4ce 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -1761,15 +1761,13 @@
                      description="Potential temperature"/>
 
                 <var name="temperature" type="real" dimensions="nVertLevels nCells Time" units="K"
-                     description="temperature"
-                     packages="jedi_da"/>
+                     description="temperature"/>
 
                 <var name="relhum" type="real" dimensions="nVertLevels nCells Time" units="percent"
                      description="Relative humidity"/>
 
                 <var name="spechum" type="real" dimensions="nVertLevels nCells Time" units="kg kg^{-1}"
-                     description="Specific humidity"
-                     packages="jedi_da"/>
+                     description="Specific humidity"/>
 
                 <var name="v" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
                      description="Horizontal tangential velocity at edges"/>
diff --git a/src/core_atmosphere/diagnostics/Makefile b/src/core_atmosphere/diagnostics/Makefile
index 614bc1c137..4cd1865e55 100644
--- a/src/core_atmosphere/diagnostics/Makefile
+++ b/src/core_atmosphere/diagnostics/Makefile
@@ -10,6 +10,7 @@ DIAGNOSTIC_MODULES = \
 	mpas_convective_diagnostics.o \
 	mpas_pv_diagnostics.o \
 	mpas_soundings.o \
+	mpas_modellevel_diagnostics.o \
 
 mpas_isobaric_diagnostics.o: mpas_atm_diagnostics_utils.o
 
@@ -21,6 +22,7 @@ mpas_pv_diagnostics.o: mpas_atm_diagnostics_utils.o
 
 mpas_soundings.o:
 
+mpas_modellevel_diagnostics.o: mpas_atm_diagnostics_utils.o
 
 ################### Generally no need to modify below here ###################
 
diff --git a/src/core_atmosphere/diagnostics/Registry_diagnostics.xml b/src/core_atmosphere/diagnostics/Registry_diagnostics.xml
index b9e7dc5682..0a9c57e7d6 100644
--- a/src/core_atmosphere/diagnostics/Registry_diagnostics.xml
+++ b/src/core_atmosphere/diagnostics/Registry_diagnostics.xml
@@ -19,6 +19,9 @@
 <!-- soundings -->
 #include "Registry_soundings.xml"
 
+<!-- additional model-level diagnostics -->
+#include "Registry_modellevel.xml"
+
 <!-- ******************************* -->
 <!-- End includes from diagnostics -->
 <!-- ******************************* -->
diff --git a/src/core_atmosphere/diagnostics/Registry_modellevel.xml b/src/core_atmosphere/diagnostics/Registry_modellevel.xml
new file mode 100644
index 0000000000..e6117845c6
--- /dev/null
+++ b/src/core_atmosphere/diagnostics/Registry_modellevel.xml
@@ -0,0 +1,15 @@
+<!-- ******************************************* -->
+<!-- Diagnostics for model-level height surfaces -->
+<!-- ******************************************* -->
+
+<!-- temperature and spechum are already added in core_atmosphere/Registry.xml -->
+
+<!-- <var_struct name="diag" time_levs="1">
+
+        <var name="temperature" type="real" dimensions="nVertLevels nCells Time" units="K"
+             description="temperature" />
+
+        <var name="spechum" type="real" dimensions="nVertLevels nCells Time" units="kg kg^{-1}"
+             description="specific humidity"/>
+
+</var_struct> -->
diff --git a/src/core_atmosphere/diagnostics/mpas_atm_diagnostics_manager.F b/src/core_atmosphere/diagnostics/mpas_atm_diagnostics_manager.F
index fb57411d1d..2f04f3a010 100644
--- a/src/core_atmosphere/diagnostics/mpas_atm_diagnostics_manager.F
+++ b/src/core_atmosphere/diagnostics/mpas_atm_diagnostics_manager.F
@@ -38,6 +38,7 @@ subroutine mpas_atm_diag_setup(stream_mgr, configs, structs, clock, dminfo)
         use mpas_convective_diagnostics, only : convective_diagnostics_setup
         use mpas_pv_diagnostics, only : pv_diagnostics_setup
         use mpas_soundings, only : soundings_setup
+        use mpas_modellevel_diagnostics, only : modellevel_diagnostics_setup
 
         implicit none
 
@@ -59,6 +60,7 @@ subroutine mpas_atm_diag_setup(stream_mgr, configs, structs, clock, dminfo)
         call convective_diagnostics_setup(structs, clock)
         call pv_diagnostics_setup(structs, clock)
         call soundings_setup(configs, structs, clock, dminfo)
+        call modellevel_diagnostics_setup(structs, clock)
 
     end subroutine mpas_atm_diag_setup
 
@@ -105,6 +107,7 @@ subroutine mpas_atm_diag_compute()
         use mpas_convective_diagnostics, only : convective_diagnostics_compute
         use mpas_pv_diagnostics, only : pv_diagnostics_compute
         use mpas_soundings, only : soundings_compute
+        use mpas_modellevel_diagnostics, only : modellevel_diagnostics_compute
 
         implicit none
 
@@ -115,6 +118,7 @@ subroutine mpas_atm_diag_compute()
         call convective_diagnostics_compute()
         call pv_diagnostics_compute()
         call soundings_compute()
+        call modellevel_diagnostics_compute()
 
     end subroutine mpas_atm_diag_compute
 
diff --git a/src/core_atmosphere/diagnostics/mpas_modellevel_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_modellevel_diagnostics.F
new file mode 100644
index 0000000000..799133fb4a
--- /dev/null
+++ b/src/core_atmosphere/diagnostics/mpas_modellevel_diagnostics.F
@@ -0,0 +1,146 @@
+! Copyright (c) 2022, University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at https://mpas-dev.github.io/license.html
+!
+module mpas_modellevel_diagnostics
+
+    use mpas_derived_types, only : MPAS_pool_type, MPAS_clock_type
+    use mpas_kind_types, only : RKIND
+    use mpas_constants, only : rvord
+
+    type (MPAS_pool_type), pointer :: mesh
+    type (MPAS_pool_type), pointer :: state
+    type (MPAS_pool_type), pointer :: diag
+    !type (MPAS_pool_type), pointer :: diag_physics
+
+    type (MPAS_clock_type), pointer :: clock
+
+    public :: modellevel_diagnostics_setup, &
+              modellevel_diagnostics_compute
+
+    private
+
+
+    contains
+
+
+    !-----------------------------------------------------------------------
+    !  routine modellevel_diagnostics_setup
+    !
+    !> \brief Initialize the modellevel diagnostic module
+    !> \author Jihyeon Jang
+    !> \date   30 January 2026
+    !> \details
+    !>  Initialize the diagnostic and save pointers to subpools for
+    !>  reuse in this module
+    !
+    !-----------------------------------------------------------------------
+    subroutine modellevel_diagnostics_setup(all_pools, simulation_clock)
+
+        use mpas_derived_types, only : MPAS_pool_type, MPAS_clock_type
+        use mpas_pool_routines, only : mpas_pool_get_subpool
+
+        implicit none
+
+        type (MPAS_pool_type), pointer :: all_pools
+        type (MPAS_clock_type), pointer :: simulation_clock
+
+        clock => simulation_clock
+
+        call mpas_pool_get_subpool(all_pools, 'mesh', mesh)
+        call mpas_pool_get_subpool(all_pools, 'state', state)
+        call mpas_pool_get_subpool(all_pools, 'diag', diag)
+
+        !call mpas_pool_get_array(diag, 'temperature', temperature)
+        !call mpas_pool_get_array(diag, 'spechum', spechum)
+        !
+        ! Zero-out the initial field
+        !
+        !temperature(:,:) = 0.0_RKIND
+        !spechum(:,:) = 0.0_RKIND
+
+    end subroutine modellevel_diagnostics_setup
+
+
+    !-----------------------------------------------------------------------
+    !  routine modellevel_diagnostics_compute
+    !
+    !> \brief Compute diagnostic before model output is written
+    !> \author Jihyeon Jang
+    !> \date   30 January 2026
+    !> \details
+    !>  Compute diagnostic before model output is written
+    !>  The following fields are computed by this routine:
+    !>     temperature
+    !>     spechum
+    !
+    !-----------------------------------------------------------------------
+    subroutine modellevel_diagnostics_compute()
+
+        use mpas_atm_diagnostics_utils, only : MPAS_field_will_be_written
+        use mpas_pool_routines, only : mpas_pool_get_dimension, mpas_pool_get_array
+
+        implicit none
+
+        integer :: iCell, k
+        integer :: time_lev
+        integer, pointer :: nCellsSolve, nVertLevels
+        integer, pointer :: index_qv
+
+        real (kind=RKIND), dimension(:,:), pointer :: temperature
+        real (kind=RKIND), dimension(:,:), pointer :: spechum
+        real (kind=RKIND), dimension(:,:), pointer :: exner, theta_m
+        real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+
+        logical :: need_ml_diags, need_temperature, need_spechum
+
+        time_lev = 1
+
+        need_ml_diags = .false.
+        need_temperature = MPAS_field_will_be_written('temperature')
+        need_ml_diags = need_ml_diags .or. need_temperature
+        need_spechum = MPAS_field_will_be_written('spechum')
+        need_ml_diags = need_ml_diags .or. need_spechum
+
+
+        if (need_ml_diags) then
+           call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve)
+           call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
+        end if
+
+        if (need_temperature .or. need_spechum) then
+           call mpas_pool_get_dimension(state, 'index_qv', index_qv)
+           call mpas_pool_get_array(state, 'scalars', scalars, time_lev)
+        end if
+
+        if (need_temperature) then
+           call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev)
+           call mpas_pool_get_array(diag, 'exner', exner)
+           call mpas_pool_get_array(diag, 'temperature', temperature)
+        end if
+
+        if (need_spechum) then
+           call mpas_pool_get_array(diag, 'spechum', spechum)
+        end if
+
+        if (need_temperature) then
+           do iCell=1,nCellsSolve
+           do k=1,nVertLevels
+              temperature(k,iCell) = (theta_m(k,iCell)/(1._RKIND+rvord*scalars(index_qv,k,iCell)))*exner(k,iCell)
+           end do
+           end do
+        end if
+
+        if (need_spechum) then
+           do iCell=1,nCellsSolve
+           do k=1,nVertLevels
+              spechum(k,iCell) = scalars(index_qv,k,iCell) / (1.0_RKIND+scalars(index_qv,k,iCell))
+           end do
+           end do
+        end if
+
+    end subroutine modellevel_diagnostics_compute
+
+end module mpas_modellevel_diagnostics