diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
index 63ddc6f0..82228032 100644
--- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
+++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
@@ -1,12 +1,17 @@
+! SPDX-FileCopyrightText: Copyright (C) 2024 U.S. National Science Foundation National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+!> This module, the MPAS subdriver, manages the life cycle (i.e., initialization, running, and
+!> finalization) of MPAS as a dynamical core within CAM-SIMA as well as potentially other
+!> host models.
+!>
+!> It is a ground-up implementation that not only adheres to the Fortran 2018 standard, but also
+!> incorporates a modern object-oriented design. As such, the implementation details of MPAS are
+!> abstracted away from CAM-SIMA, which enables a more stable interface between the two.
+!>
+!> Users should begin by creating an "instance" of MPAS dynamical core from the `mpas_dynamical_core_type`
+!> derived type. Then, interaction with the instance is done through its public type-bound procedures.
 module dyn_mpas_subdriver
-    !-------------------------------------------------------------------------------
-    ! module dyn_mpas_subdriver
-    !
-    ! This module manages the life cycle (i.e., initialization, running, and
-    ! finalization) of MPAS as a dynamical core within CAM-SIMA.
-    !
-    !-------------------------------------------------------------------------------
-
     use, intrinsic :: iso_fortran_env, only: output_unit
     ! Module(s) from external libraries.
 #ifdef MPAS_USE_MPI_F08
@@ -26,7 +31,8 @@ module dyn_mpas_subdriver
     public :: mpas_dynamical_core_type
 
     abstract interface
-        ! This interface is compatible with `endrun` from CAM-SIMA.
+        !> This procedure interface is modeled after the `endrun` subroutine from CAM-SIMA.
+        !> It will be called whenever MPAS dynamical core encounters a fatal error and cannot continue.
         subroutine model_error_if(message, file, line)
             character(*),           intent(in) :: message
             character(*), optional, intent(in) :: file
@@ -51,8 +57,9 @@ end subroutine model_error_if
     integer, parameter :: mpas_dynamical_core_real_kind = rkind
 
     !> The "class" of MPAS dynamical core.
-    !> Important data structures like states of MPAS dynamical core are encapsulated inside this derived type to prevent misuse.
-    !> Type-bound procedures provide well-defined APIs for CAM-SIMA to interact with MPAS dynamical core.
+    !> Important data structures like the internal states of MPAS dynamical core are encapsulated inside this derived type
+    !> to prevent misuse. Type-bound procedures provide stable and well-defined APIs for CAM-SIMA to interact with
+    !> MPAS dynamical core.
     type :: mpas_dynamical_core_type
         private
 
@@ -67,7 +74,7 @@ end subroutine model_error_if
         integer :: mpi_rank = 0
         logical :: mpi_rank_root = .false.
 
-        ! Actual implementation is supplied at runtime.
+        ! Actual implementation is supplied at run-time.
         procedure(model_error_if), nopass, pointer :: model_error => null()
 
         type(core_type), pointer :: corelist => null()
@@ -104,7 +111,7 @@ end subroutine model_error_if
         procedure, pass, public :: run => dyn_mpas_run
         procedure, pass, public :: final => dyn_mpas_final
 
-        ! Accessor subroutines for users to access internal states of MPAS dynamical core.
+        ! Accessor procedures for users to access the internal states of MPAS dynamical core.
 
         procedure, pass, public :: get_constituent_name => dyn_mpas_get_constituent_name
         procedure, pass, public :: get_constituent_index => dyn_mpas_get_constituent_index
@@ -161,9 +168,13 @@ end subroutine model_error_if
 
     !> This derived type conveys information similar to the `var` and `var_array` elements in MPAS registry.
     !> For example, in MPAS registry, the "xCell" variable is described as:
-    !>     <var name="xCell" type="real" dimensions="nCells" units="m" description="Cartesian x-coordinate of cells" />
+    !> ```
+    !> <var name="xCell" type="real" dimensions="nCells" units="m" description="Cartesian x-coordinate of cells" />
+    !> ```
     !> Here, it is described as:
-    !>     var_info_type(name="xCell", type="real", rank=1)
+    !> ```
+    !> var_info_type(name="xCell", type="real", rank=1)
+    !> ```
     !> However, note that MPAS treats the "Time" dimension specially. It is implemented as 1-d pointer arrays of
     !> custom derived types. For a variable with the "Time" dimension, its rank needs to be subtracted by one.
     type :: var_info_type
@@ -335,7 +346,7 @@ end subroutine dyn_mpas_debug_print
     !> If `value` contains more than one element, the elements will be stringified, delimited by `separator`, then concatenated.
     !> If `value` contains exactly one element, the element will be stringified without using `separator`.
     !> If `value` contains zero element or is of unsupported data types, an empty character string is produced.
-    !> If `separator` is not supplied, it defaults to `, ` (i.e., a comma and a space).
+    !> If `separator` is not supplied, it defaults to ", " (i.e., a comma and a space).
     !> (KCW, 2024-02-04)
     pure function stringify(value, separator)
         use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
@@ -435,13 +446,12 @@ end function stringify
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_phase1
     !
-    !> \brief  Tracks `mpas_init` up to the point of reading namelist
+    !> \brief  Track `mpas_init` up to the point of reading namelist.
     !> \author Michael Duda
     !> \date   19 April 2019
     !> \details
     !>  This subroutine follows the stand-alone MPAS subdriver up to, but not
     !>  including, the point where namelist is read.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-02-02)
     !
     !-------------------------------------------------------------------------------
@@ -530,9 +540,9 @@ subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_level, log
         ! initialization steps.
         !
         ! We need:
-        ! 1) `domain_ptr` to be allocated;
-        ! 2) `dmpar_init` to be completed for accessing `dminfo`;
-        ! 3) `*_setup_core` to assign the `setup_log` procedure pointer.
+        ! 1. `domain_ptr` to be allocated;
+        ! 2. `dmpar_init` to be completed for accessing `dminfo`;
+        ! 3. `*_setup_core` to assign the `setup_log` procedure pointer.
         ierr = self % domain_ptr % core % setup_log(self % domain_ptr % loginfo, self % domain_ptr, unitnumbers=mpas_log_unit)
 
         if (ierr /= 0) then
@@ -540,20 +550,20 @@ subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_level, log
                 subname, __LINE__)
         end if
 
-        ! At this point, we should be ready to read namelist in `dyn_comp::dyn_readnl`.
+        ! At this point, we should be ready to read namelist in `dyn_mpas_read_namelist`.
         call self % debug_print(log_level_debug, subname // ' completed')
     end subroutine dyn_mpas_init_phase1
 
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_read_namelist
     !
-    !> \brief  Tracks `mpas_init` where namelist is being read
+    !> \brief  Track `mpas_init` where namelist is being read.
     !> \author Kuan-Chih Wang
     !> \date   2024-02-09
     !> \details
     !>  This subroutine calls upstream MPAS functionality for reading its own
     !>  namelist. After that, override designated namelist variables according to
-    !>  information provided from CAM-SIMA.
+    !>  the information provided from CAM-SIMA.
     !
     !-------------------------------------------------------------------------------
     subroutine dyn_mpas_read_namelist(self, namelist_path, &
@@ -591,13 +601,13 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
                 subname, __LINE__)
         end if
 
-        ! Override designated namelist variables according to information provided from CAM-SIMA.
-        ! These include runtime settings that cannot be determined beforehand.
+        ! Override designated namelist variables according to the information provided from CAM-SIMA.
+        ! These include run-time settings that cannot be determined beforehand.
 
         call self % debug_print(log_level_info, 'Overriding designated namelist variables')
 
         ! CAM-SIMA seems to follow "NetCDF Climate and Forecast (CF) Metadata Conventions" for calendar names. See
-        ! CF-1.11, section "4.4.1. Calendar".
+        ! CF-1.12, section "4.4.2. Calendar".
         ! However, this is not the case for MPAS. Translate calendar names between CF and MPAS.
         select case (trim(adjustl(cf_calendar)))
             case ('360_day')
@@ -605,7 +615,7 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
             case ('365_day', 'noleap')
                 mpas_calendar = 'gregorian_noleap'
             case ('gregorian', 'standard')
-                ! `gregorian` is a deprecated alternative name for `standard`.
+                ! "gregorian" is a deprecated alternative name for "standard".
                 mpas_calendar = 'gregorian'
             case default
                 call self % model_error('Unsupported calendar type "' // trim(adjustl(cf_calendar)) // '"', &
@@ -618,9 +628,9 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
         call self % debug_print(log_level_debug, 'config_calendar_type = ' // stringify([config_pointer_c]))
         nullify(config_pointer_c)
 
-        ! MPAS represents date and time in ISO 8601 format. However, the separator between date and time is `_`
-        ! instead of standard `T`.
-        ! Format in `YYYY-MM-DD_hh:mm:ss` is acceptable.
+        ! MPAS represents date and time in ISO 8601 format. However, the separator between date and time is "_"
+        ! instead of standard "T".
+        ! Format in "YYYY-MM-DD_hh:mm:ss" is acceptable.
         call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_start_time')
 
         config_pointer_c = stringify(start_date_time(1:3), '-') // '_' // stringify(start_date_time(4:6), ':')
@@ -633,7 +643,7 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
         call self % debug_print(log_level_debug, 'config_stop_time = ' // stringify([config_pointer_c]))
         nullify(config_pointer_c)
 
-        ! Format in `DD_hh:mm:ss` is acceptable.
+        ! Format in "DD_hh:mm:ss" is acceptable.
         call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_run_duration')
 
         config_pointer_c = stringify([run_duration(1)]) // '_' // stringify(run_duration(2:4), ':')
@@ -654,20 +664,21 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
         call self % debug_print(log_level_debug, 'config_do_restart = ' // stringify([config_pointer_l]))
         nullify(config_pointer_l)
 
+        ! At this point, we should be ready to follow up with the rest of MPAS framework initialization
+        ! in `dyn_mpas_init_phase2`.
         call self % debug_print(log_level_debug, subname // ' completed')
     end subroutine dyn_mpas_read_namelist
 
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_phase2
     !
-    !> \brief  Tracks `mpas_init` after namelist has been read
+    !> \brief  Track `mpas_init` after namelist has been read.
     !> \author Michael Duda
     !> \date   19 April 2019
     !> \details
     !>  This subroutine follows the stand-alone MPAS subdriver from the point
     !>  where we call the second phase of MPAS framework initialization up
-    !>  to the check on the existence of the `streams.<core>` file.
-    !> \addenda
+    !>  to the check on the existence of the "streams.<core>" file.
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-02-07)
     !
     !-------------------------------------------------------------------------------
@@ -704,7 +715,7 @@ subroutine dyn_mpas_init_phase2(self, pio_iosystem)
         ! Initialize MPAS framework with the supplied PIO system descriptor.
         call mpas_framework_init_phase2(self % domain_ptr, io_system=pio_iosystem)
 
-        ! Instantiate `streaminfo` but do not actually initialize it. Any queries made to it will always return `.false.`.
+        ! Instantiate `streaminfo`, but do not actually initialize it. Any queries made to it will always return `.false.`.
         ! This is the intended behavior because MPAS as a dynamical core is not responsible for managing IO.
         self % domain_ptr % streaminfo => mpas_stream_inquiry_new_streaminfo()
 
@@ -752,25 +763,24 @@ subroutine dyn_mpas_init_phase2(self, pio_iosystem)
         end if
 
         ! At this point, we should be ready to set up decompositions, build halos, allocate blocks, etc.
-        ! in `dyn_grid::model_grid_init`.
+        ! in `dyn_mpas_init_phase3`.
         call self % debug_print(log_level_debug, subname // ' completed')
     end subroutine dyn_mpas_init_phase2
 
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_phase3
     !
-    !> \brief  Tracks `mpas_init` up to the point of calling `atm_core_init`
+    !> \brief  Track `mpas_init` up to the point of calling `atm_core_init`.
     !> \author Michael Duda
     !> \date   19 April 2019
     !> \details
     !>  This subroutine follows the stand-alone MPAS subdriver after the check on
-    !>  the existence of the `streams.<core>` file up to, but not including,
+    !>  the existence of the "streams.<core>" file up to, but not including,
     !>  the point where `atm_core_init` is called. It completes MPAS framework
     !>  initialization, including the allocation of all blocks and fields managed
-    !>  by MPAS. However, scalars are allocated but not yet defined.
+    !>  by MPAS. However, note that scalars are allocated, but not yet defined.
     !>  `dyn_mpas_define_scalar` must be called afterwards. Also note that MPAS uses
     !>  the term "scalar", but CAM-SIMA calls it "constituent".
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-03-06)
     !
     !-------------------------------------------------------------------------------
@@ -804,7 +814,7 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file)
 
         call self % debug_print(log_level_info, 'Number of constituents is ' // stringify([self % number_of_constituents]))
 
-        ! Adding a config named `cam_pcnst` with the number of constituents will indicate to MPAS that
+        ! Adding a config named "cam_pcnst" with the number of constituents will indicate to MPAS that
         ! it is operating as a dynamical core, and therefore it needs to allocate scalars separately
         ! from other Registry-defined fields. The special logic is located in `atm_setup_block`.
         ! This must be done before calling `mpas_bootstrap_framework_phase1`.
@@ -834,8 +844,8 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file)
         ! Finish setting up fields.
         call mpas_bootstrap_framework_phase2(self % domain_ptr, pio_file_desc=pio_file)
 
-        ! `num_scalars` is a dimension variable, but it only exists in MPAS `state` pool.
-        ! Fix this inconsistency by also adding it to MPAS `dimension` pool.
+        ! "num_scalars" is a dimension variable, but it only exists in MPAS "state" pool.
+        ! Fix this inconsistency by also adding it to MPAS "dimension" pool.
         call self % get_pool_pointer(mpas_pool, 'state')
 
         call mpas_pool_get_dimension(mpas_pool, 'num_scalars', num_scalars)
@@ -854,13 +864,22 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file)
         nullify(mpas_pool)
         nullify(num_scalars)
 
+        ! At this point, what follows next depends on the specific use case. In no particular order:
+        ! * Use `dyn_mpas_define_scalar` to define the names of constituents at run-time.
+        ! * Use `dyn_mpas_read_write_stream` to read mesh variables.
+        !   * Follow up with a call to `dyn_mpas_compute_unit_vector` immediately. This is by design.
+        ! * For setting analytic initial condition, use `get_variable_pointer` to inject data directly into MPAS memory.
+        !   * Use `dyn_mpas_compute_edge_wind` where appropriate.
+        !   * Use `dyn_mpas_exchange_halo` where appropriate.
+        ! * Use `dyn_mpas_read_write_stream` to read initial condition or restart.
+        ! * Finally, use `dyn_mpas_init_phase4` to conclude the initialization of MPAS dynamical core.
         call self % debug_print(log_level_debug, subname // ' completed')
     end subroutine dyn_mpas_init_phase3
 
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_define_scalar
     !
-    !> \brief  Defines the names of constituents at run-time
+    !> \brief  Define the names of constituents at run-time.
     !> \author Michael Duda
     !> \date   21 May 2020
     !> \details
@@ -873,7 +892,6 @@ end subroutine dyn_mpas_init_phase3
     !>  this constrain. Index mapping between MPAS scalars and constituent names
     !>  can be looked up through `index_constituent_to_mpas_scalar` and
     !>  `index_mpas_scalar_to_constituent`.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-05-19)
     !
     !-------------------------------------------------------------------------------
@@ -886,9 +904,9 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
         character(*), intent(in) :: constituent_name(:)
         logical, intent(in) :: is_water_species(:)
 
-        ! Possible CCPP standard names of `qv`, which denotes water vapor mixing ratio.
-        ! They are hard-coded here because MPAS needs to know where `qv` is.
-        ! Index 1 is exactly what MPAS wants. Others also work, but need to be converted.
+        !> Possible CCPP standard names of `qv`, which denotes water vapor mixing ratio.
+        !> They are hard-coded here because MPAS needs to know where `qv` is.
+        !> Index 1 is exactly what MPAS wants. Others also work, but need to be converted.
         character(*), parameter :: mpas_scalar_qv_standard_name(*) = [ character(strkind) :: &
             'water_vapor_mixing_ratio_wrt_dry_air', &
             'water_vapor_mixing_ratio_wrt_moist_air', &
@@ -981,7 +999,7 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
 
             if (index_qv == 0) then
                 call self % model_error('Constituent names must contain one of: ' // &
-                    stringify(mpas_scalar_qv_standard_name) // ' and it must be a water species', subname, __LINE__)
+                    stringify(mpas_scalar_qv_standard_name) // ', and it must be a water species', subname, __LINE__)
             end if
         end if
 
@@ -1059,7 +1077,7 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
         call mpas_pool_add_dimension(mpas_pool, 'moist_start', index_water_start)
         call mpas_pool_add_dimension(mpas_pool, 'moist_end', index_water_end)
 
-        ! MPAS `state` pool has two time levels.
+        ! MPAS "state" pool has two time levels.
         time_level = 2
 
         do i = 1, time_level
@@ -1100,7 +1118,7 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
         call mpas_pool_add_dimension(mpas_pool, 'moist_start', index_water_start)
         call mpas_pool_add_dimension(mpas_pool, 'moist_end', index_water_end)
 
-        ! MPAS `tend` pool only has one time level.
+        ! MPAS "tend" pool only has one time level.
         time_level = 1
 
         do i = 1, time_level
@@ -1131,7 +1149,7 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
 
         nullify(mpas_pool)
 
-        ! For consistency, also add dimension variables to MPAS `dimension` pool.
+        ! For consistency, also add dimension variables to MPAS "dimension" pool.
 
         call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'index_qv', index_qv)
         call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_start', index_water_start)
@@ -1147,7 +1165,7 @@ end subroutine dyn_mpas_define_scalar
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_read_write_stream
     !
-    !> \brief  Read or write an MPAS stream
+    !> \brief  Read or write an MPAS stream.
     !> \author Kuan-Chih Wang
     !> \date   2024-03-15
     !> \details
@@ -1259,7 +1277,7 @@ end subroutine dyn_mpas_read_write_stream
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_stream_with_pool
     !
-    !> \brief  Initialize an MPAS stream with an accompanying MPAS pool
+    !> \brief  Initialize an MPAS stream with an accompanying MPAS pool.
     !> \author Kuan-Chih Wang
     !> \date   2024-03-14
     !> \details
@@ -1299,9 +1317,9 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_stream_with_pool'
         character(strkind) :: stream_filename
         integer :: i, ierr, stream_format
-        ! Whether a variable is present on the file (i.e., `pio_file`).
+        !> Whether a variable is present on the file (i.e., `pio_file`).
         logical, allocatable :: var_is_present(:)
-        ! Whether a variable is type, kind and rank compatible with what MPAS expects on the file (i.e., `pio_file`).
+        !> Whether a variable is type, kind, and rank compatible with what MPAS expects on the file (i.e., `pio_file`).
         logical, allocatable :: var_is_tkr_compatible(:)
         type(field0dchar), pointer :: field_0d_char
         type(field1dchar), pointer :: field_1d_char
@@ -1637,11 +1655,11 @@ subroutine add_stream_attribute_0d(attribute_name, attribute_value)
                         trim(adjustl(attribute_name)), attribute_value, syncval=.false., ierr=ierr)
                 type is (logical)
                     if (attribute_value) then
-                        ! Logical `.true.` becomes character string `YES`.
+                        ! Logical `.true.` becomes character string "YES".
                         call mpas_writestreamatt(mpas_stream, &
                             trim(adjustl(attribute_name)), 'YES', syncval=.false., ierr=ierr)
                     else
-                        ! Logical `.false.` becomes character string `NO`.
+                        ! Logical `.false.` becomes character string "NO".
                         call mpas_writestreamatt(mpas_stream, &
                             trim(adjustl(attribute_name)), 'NO', syncval=.false., ierr=ierr)
                     end if
@@ -1691,10 +1709,12 @@ end subroutine add_stream_attribute_1d
     end subroutine dyn_mpas_init_stream_with_pool
 
     !> Parse a stream name, which consists of one or more stream name fragments, and return the corresponding variable information
-    !> as a list of `var_info_type`. Multiple stream name fragments should be separated by `+` (i.e., a plus, meaning "addition"
-    !> operation) or `-` (i.e., a minus, meaning "subtraction" operation).
-    !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", "restart") or a single variable name.
-    !> Duplicate variable names in the resulting list are discarded.
+    !> as a list of `var_info_type`. Multiple stream name fragments should be separated by "+" (i.e., a plus, meaning "addition"
+    !> operation) or "-" (i.e., a minus, meaning "subtraction" operation).
+    !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", etc.) or a single variable name.
+    !> For example, a stream name of "invariant+input+restart" means the union of variables in the "invariant", "input", and
+    !> "restart" streams.
+    !> Duplicate variable information in the resulting list is discarded.
     !> (KCW, 2024-06-01)
     pure function parse_stream_name(stream_name) result(var_info_list)
         character(*), intent(in) :: stream_name
@@ -1782,7 +1802,7 @@ pure function parse_stream_name(stream_name) result(var_info_list)
     end function parse_stream_name
 
     !> Parse a stream name fragment and return the corresponding variable information as a list of `var_info_type`.
-    !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", "restart") or a single variable name.
+    !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", etc.) or a single variable name.
     !> (KCW, 2024-06-01)
     pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_list)
         character(*), intent(in) :: stream_name_fragment
@@ -1837,6 +1857,7 @@ end function parse_stream_name_fragment
 
     !> Return the index of unique elements in `array`, which can be any intrinsic data types, as an integer array.
     !> If `array` contains zero element or is of unsupported data types, an empty integer array is produced.
+    !> For example, `index_unique([1, 2, 3, 1, 2, 3, 4, 5])` returns `[1, 2, 3, 7, 8]`.
     !> (KCW, 2024-03-22)
     pure function index_unique(array)
         use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
@@ -1916,13 +1937,13 @@ end function index_unique
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_check_variable_status
     !
-    !> \brief  Check and return variable status on the given file
+    !> \brief  Check and return variable status on the given file.
     !> \author Kuan-Chih Wang
     !> \date   2024-06-04
     !> \details
     !>  On the given file (i.e., `pio_file`), this subroutine checks whether the
     !>  given variable (i.e., `var_info`) is present, and whether it is "TKR"
-    !>  compatible with what MPAS expects. "TKR" means type, kind and rank.
+    !>  compatible with what MPAS expects. "TKR" means type, kind, and rank.
     !>  This subroutine can handle both ordinary variables and variable arrays.
     !>  They are indicated by the `var` and `var_array` elements, respectively,
     !>  in MPAS registry. For an ordinary variable, the checks are performed on
@@ -2348,13 +2369,12 @@ end subroutine dyn_mpas_check_variable_status
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_exchange_halo
     !
-    !> \brief  Updates the halo layers of the named field
+    !> \brief  Update the halo layers of the named field.
     !> \author Michael Duda
     !> \date   16 January 2020
     !> \details
     !>  Given a field name that is defined in MPAS registry, this subroutine updates
     !>  the halo layers for that field.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-03-18)
     !
     !-------------------------------------------------------------------------------
@@ -2531,21 +2551,21 @@ end subroutine dyn_mpas_exchange_halo
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_compute_unit_vector
     !
-    !> \brief  Computes local east, north and edge-normal unit vectors
+    !> \brief  Compute local east, north, and edge-normal unit vectors.
     !> \author Michael Duda
     !> \date   15 January 2020
     !> \details
     !>  This subroutine computes the local east and north unit vectors at all cells,
-    !>  storing the results in MPAS `mesh` pool as `east` and `north`, respectively.
-    !>  It also computes the edge-normal unit vectors at all edges by calling
-    !>  `mpas_initialize_vectors`. Before calling this subroutine, MPAS `mesh` pool
-    !>  must contain `latCell` and `lonCell` that are valid for all cells (not just
-    !>  solve cells), plus any additional variables that are required by
+    !>  storing the results in MPAS "mesh" pool as the "east" and "north" variables,
+    !>  respectively. It also computes the edge-normal unit vectors at all edges by
+    !>  calling `mpas_initialize_vectors`.
+    !>  Before calling this subroutine, MPAS "mesh" pool must contain the "latCell"
+    !>  and "lonCell" variables that are valid for all cells (not just solve cells),
+    !>  plus any additional variables that are required by
     !>  `mpas_initialize_vectors`.
     !>  For stand-alone MPAS, the whole deal is handled by `init_dirs_forphys`
     !>  during physics initialization. However, MPAS as a dynamical core does
-    !>  not have physics, hence this subroutine.
-    !> \addenda
+    !>  not have physics, hence the existence of this subroutine.
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-04-23)
     !
     !-------------------------------------------------------------------------------
@@ -2613,20 +2633,20 @@ end subroutine dyn_mpas_compute_unit_vector
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_compute_edge_wind
     !
-    !> \brief  Computes the edge-normal wind vectors at edge points
+    !> \brief  Compute the edge-normal wind (tendency) vectors at edge points.
     !> \author Michael Duda
     !> \date   16 January 2020
     !> \details
-    !>  This subroutine computes the edge-normal wind vectors at edge points
-    !>  (i.e., `u` in MPAS `state` pool) from the wind components at cell points
-    !>  (i.e., `uReconstruct{Zonal,Meridional}` in MPAS `diag` pool). In MPAS, the
-    !>  former are PROGNOSTIC variables, while the latter are DIAGNOSTIC variables
-    !>  that are "reconstructed" from the former. This subroutine is essentially the
-    !>  inverse function of that reconstruction. The purpose is to provide an
-    !>  alternative way for MPAS to initialize from zonal and meridional wind
-    !>  components at cell points. If `wind_tendency` is `.true.`, this subroutine
-    !>  operates on the wind tendency due to physics instead.
-    !> \addenda
+    !>  This subroutine computes the edge-normal wind vectors at edge points (i.e.,
+    !>  the "u" variable in MPAS "state" pool) from the wind components at cell
+    !>  points (i.e., the "uReconstruct{Zonal,Meridional}" variables in MPAS "diag"
+    !>  pool). In MPAS, the former are PROGNOSTIC variables, while the latter are
+    !>  DIAGNOSTIC variables that are "reconstructed" from the former.
+    !>  This subroutine is essentially the inverse function of that reconstruction.
+    !>  The purpose is to provide an alternative way for MPAS to initialize from
+    !>  zonal and meridional wind components at cell points.
+    !>  If `wind_tendency` is `.true.`, this subroutine operates on the wind
+    !>  tendency due to physics instead.
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-05-08)
     !
     !-------------------------------------------------------------------------------
@@ -2734,14 +2754,13 @@ end subroutine dyn_mpas_compute_edge_wind
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_phase4
     !
-    !> \brief  Tracks `atm_core_init` to finish MPAS dynamical core initialization
+    !> \brief  Track `atm_core_init` to finish MPAS dynamical core initialization.
     !> \author Michael Duda
     !> \date   29 February 2020
     !> \details
     !>  This subroutine completes MPAS dynamical core initialization.
     !>  Essentially, it closely follows what is done in `atm_core_init`, but without
     !>  any calls to MPAS diagnostics manager or MPAS stream manager.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-05-25)
     !
     !-------------------------------------------------------------------------------
@@ -2857,7 +2876,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
             call self % model_error('Failed to build halo exchange groups', subname, __LINE__)
         end if
 
-        ! Variables in MPAS `state` pool have more than one time level. Copy the values from the first time level of
+        ! Variables in MPAS "state" pool have more than one time level. Copy the values from the first time level of
         ! such variables into all subsequent time levels to initialize them.
         call self % get_variable_pointer(config_do_restart, 'cfg', 'config_do_restart')
 
@@ -3018,7 +3037,7 @@ end subroutine dyn_mpas_init_phase4
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_run
     !
-    !> \brief  Integrates the dynamical states with time
+    !> \brief  Integrate the dynamical states with time.
     !> \author Michael Duda
     !> \date   29 February 2020
     !> \details
@@ -3027,7 +3046,6 @@ end subroutine dyn_mpas_init_phase4
     !>  the coupling time interval is reached.
     !>  Essentially, it closely follows what is done in `atm_core_run`, but without
     !>  any calls to MPAS diagnostics manager or MPAS stream manager.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-06-21)
     !
     !-------------------------------------------------------------------------------
@@ -3094,7 +3112,7 @@ subroutine dyn_mpas_run(self)
             ! Current states are in time level 1. Upon exit, time level 2 will contain updated states.
             call atm_do_timestep(self % domain_ptr, config_dt, self % number_of_time_steps)
 
-            ! MPAS `state` pool has two time levels.
+            ! MPAS "state" pool has two time levels.
             ! Swap them after advancing a time step.
             call mpas_pool_shift_time_levels(mpas_pool_state)
 
@@ -3121,7 +3139,7 @@ subroutine dyn_mpas_run(self)
 
         call self % debug_print(log_level_info, 'Time integration of MPAS dynamical core ends at ' // trim(adjustl(date_time)))
 
-        ! Compute diagnostic variables like `pressure`, `rho` and `theta` from time level 1 of MPAS `state` pool
+        ! Compute diagnostic variables like "pressure", "rho", and "theta" from time level 1 of MPAS "state" pool
         ! by calling upstream MPAS functionality.
         call atm_compute_output_diagnostics(mpas_pool_state, 1, mpas_pool_diag, mpas_pool_mesh)
 
@@ -3134,7 +3152,7 @@ end subroutine dyn_mpas_run
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_final
     !
-    !> \brief  Finalizes MPAS dynamical core as well as its framework
+    !> \brief  Finalize MPAS dynamical core as well as its framework.
     !> \author Michael Duda
     !> \date   29 February 2020
     !> \details
@@ -3144,7 +3162,6 @@ end subroutine dyn_mpas_run
     !>  Essentially, it closely follows what is done in `atm_core_finalize` and
     !>  `mpas_finalize`, except that here, there is no need to call MPAS diagnostics
     !>  manager or MPAS stream manager.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-10-10)
     !
     !-------------------------------------------------------------------------------
@@ -3282,7 +3299,7 @@ end subroutine dyn_mpas_final
     !-------------------------------------------------------------------------------
     ! function dyn_mpas_get_constituent_name
     !
-    !> \brief  Query constituent name by its index
+    !> \brief  Query constituent name by its index.
     !> \author Kuan-Chih Wang
     !> \date   2024-05-16
     !> \details
@@ -3316,7 +3333,7 @@ end function dyn_mpas_get_constituent_name
     !-------------------------------------------------------------------------------
     ! function dyn_mpas_get_constituent_index
     !
-    !> \brief  Query constituent index by its name
+    !> \brief  Query constituent index by its name.
     !> \author Kuan-Chih Wang
     !> \date   2024-05-16
     !> \details
@@ -3352,7 +3369,7 @@ end function dyn_mpas_get_constituent_index
     !-------------------------------------------------------------------------------
     ! function dyn_mpas_map_mpas_scalar_index
     !
-    !> \brief  Map MPAS scalar index from constituent index
+    !> \brief  Map MPAS scalar index from constituent index.
     !> \author Kuan-Chih Wang
     !> \date   2024-05-16
     !> \details
@@ -3386,7 +3403,7 @@ end function dyn_mpas_map_mpas_scalar_index
     !-------------------------------------------------------------------------------
     ! function dyn_mpas_map_constituent_index
     !
-    !> \brief  Map constituent index from MPAS scalar index
+    !> \brief  Map constituent index from MPAS scalar index.
     !> \author Kuan-Chih Wang
     !> \date   2024-05-16
     !> \details
@@ -3420,12 +3437,12 @@ end function dyn_mpas_map_constituent_index
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_get_local_mesh_dimension
     !
-    !> \brief  Returns local mesh dimensions
+    !> \brief  Return local mesh dimensions.
     !> \author Kuan-Chih Wang
     !> \date   2024-05-09
     !> \details
     !>  This subroutine returns local mesh dimensions, including:
-    !>  * Numbers of local mesh cells, edges, vertices and vertical levels
+    !>  * Numbers of local mesh cells, edges, vertices, and vertical levels
     !>    on each individual task, both with/without halo points.
     !
     !-------------------------------------------------------------------------------
@@ -3488,16 +3505,15 @@ end subroutine dyn_mpas_get_local_mesh_dimension
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_get_global_mesh_dimension
     !
-    !> \brief  Returns global mesh dimensions
+    !> \brief  Return global mesh dimensions.
     !> \author Michael Duda
     !> \date   22 August 2019
     !> \details
     !>  This subroutine returns global mesh dimensions, including:
-    !>  * Numbers of global mesh cells, edges, vertices and vertical levels
+    !>  * Numbers of global mesh cells, edges, vertices, and vertical levels
     !>    across all tasks.
     !>  * Maximum numbers of mesh cells and edges/vertices among all tasks.
     !>  * Sphere radius.
-    !> \addenda
     !>  Ported and refactored for CAM-SIMA. (KCW, 2024-03-25)
     !
     !-------------------------------------------------------------------------------
@@ -3593,7 +3609,7 @@ end subroutine dyn_mpas_get_pool_pointer
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_get_variable_pointer_*
     !
-    !> \brief  A family of accessor subroutines for MPAS dynamical core instance
+    !> \brief  A family of accessor subroutines for MPAS dynamical core instance.
     !> \author Kuan-Chih Wang
     !> \date   2024-03-21
     !> \details
@@ -3602,7 +3618,7 @@ end subroutine dyn_mpas_get_pool_pointer
     !>  MPAS dynamical core instance. The `get_variable_pointer` generic interface
     !>  should be used instead of the specific ones.
     !>  WARNING:
-    !>  USE OF THIS SUBROUTINE FAMILY IS HIGHLY DISCOURAGED BECAUSE INTERNAL
+    !>  USE OF THIS SUBROUTINE FAMILY IS HIGHLY DISCOURAGED BECAUSE THE INTERNAL
     !>  STATES OF MPAS DYNAMICAL CORE INSTANCE COULD BE MODIFIED THROUGH THE
     !>  RETURNED POINTER. THESE ARE UNCHARTED WATERS SO BE SURE WHAT YOU ARE
     !>  DOING.
@@ -3980,7 +3996,7 @@ end subroutine dyn_mpas_get_variable_pointer_r5
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_get_variable_value_*
     !
-    !> \brief  A family of accessor subroutines for MPAS dynamical core instance
+    !> \brief  A family of accessor subroutines for MPAS dynamical core instance.
     !> \author Kuan-Chih Wang
     !> \date   2024-03-21
     !> \details
diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90
index 5e686f25..803f281b 100644
--- a/src/dynamics/mpas/dyn_comp.F90
+++ b/src/dynamics/mpas/dyn_comp.F90
@@ -1,3 +1,13 @@
+! SPDX-FileCopyrightText: Copyright (C) 2024 U.S. National Science Foundation National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+!> This module, part of the MPAS interface, integrates MPAS dynamical core with CAM-SIMA by
+!> implementing the necessary APIs and managing their interaction.
+!>
+!> It contains the instance of MPAS dynamical core, which is used extensively throughout CAM-SIMA.
+!> It provides core functionalities such as the initialization, running, and finalization of MPAS
+!> dynamical core. Various utility procedures for debug printing, exchanging constituent states,
+!> inquiring mesh dimensions, etc. are also provided here.
 module dyn_comp
     ! Module(s) from CESM Share.
     use shr_kind_mod, only: kind_r8 => shr_kind_r8, &
@@ -43,6 +53,7 @@ module dyn_comp
     type(mpas_dynamical_core_type) :: mpas_dynamical_core
 
     ! Local and global mesh dimensions of MPAS dynamical core.
+    ! Protected module variables that can only be initialized by `dyn_inquire_mesh_dimensions`.
     integer, protected :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
     integer, protected :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max
     real(kind_dyn_mpas), protected :: sphere_radius
@@ -113,7 +124,7 @@ subroutine dyn_readnl(namelist_path)
 
         nullify(pio_iosystem)
 
-        ! Get free units for MPAS so it can write its own log files, e.g., `log.atmosphere.0000.{out,err}`.
+        ! Get free units for MPAS so it can write its own log files, e.g., "log.atmosphere.0000.{out,err}".
         log_unit(1) = shr_file_getunit()
         log_unit(2) = shr_file_getunit()
 
@@ -136,7 +147,7 @@ subroutine dyn_readnl(namelist_path)
 
         call dyn_debug_print(debugout_info, 'Reading namelist')
 
-        ! Read MPAS-related namelist variables from `namelist_path`, e.g., `atm_in`.
+        ! Read MPAS-related namelist variables from `namelist_path`, e.g., "atm_in".
         call mpas_dynamical_core % read_namelist(namelist_path, &
             cam_calendar, start_date_time, stop_date_time, run_duration, initial_run)
 
@@ -331,7 +342,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
     end subroutine dyn_init
 
     !> Check for consistency in topography data. The presence of topography file is inferred from the `pio_file` pointer.
-    !> If topography file is used, check that the `PHIS` variable, which denotes surface geopotential,
+    !> If topography file is used, check that the "PHIS" variable, which denotes surface geopotential,
     !> is consistent with the surface geometric height in MPAS.
     !> Otherwise, if topography file is not used, check that the surface geometric height in MPAS is zero.
     !> (KCW, 2024-05-10)
@@ -352,7 +363,7 @@ subroutine check_topography_data(pio_file)
         real(kind_r8), parameter :: error_tolerance = 1.0E-3_kind_r8 ! Error tolerance for consistency check.
         real(kind_r8), allocatable :: surface_geometric_height(:)    ! Computed from topography file.
         real(kind_r8), allocatable :: surface_geopotential(:)        ! Read from topography file.
-        real(kind_dyn_mpas), pointer :: zgrid(:, :)                  ! From MPAS. Geometric height (meters) at layer interfaces.
+        real(kind_dyn_mpas), pointer :: zgrid(:, :)                  ! From MPAS. Geometric height (m) at layer interfaces.
 
         call dyn_debug_print(debugout_debug, subname // ' entered')
 
@@ -416,9 +427,9 @@ subroutine set_analytic_initial_condition()
         integer, allocatable :: global_grid_index(:)
         real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :)
         real(kind_r8), allocatable :: lat_rad(:), lon_rad(:)
-        real(kind_r8), allocatable :: z_int(:, :)   ! Geometric height (meters) at layer interfaces.
+        real(kind_r8), allocatable :: z_int(:, :)   ! Geometric height (m) at layer interfaces.
                                                     ! Dimension and vertical index orders follow CAM-SIMA convention.
-        real(kind_dyn_mpas), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces.
+        real(kind_dyn_mpas), pointer :: zgrid(:, :) ! Geometric height (m) at layer interfaces.
                                                     ! Dimension and vertical index orders follow MPAS convention.
 
         call dyn_debug_print(debugout_debug, subname // ' entered')
diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90
index eeca26e2..2f5cd7ac 100644
--- a/src/dynamics/mpas/dyn_coupling.F90
+++ b/src/dynamics/mpas/dyn_coupling.F90
@@ -1,3 +1,10 @@
+! SPDX-FileCopyrightText: Copyright (C) 2024 U.S. National Science Foundation National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+!> This module, part of the MPAS interface, integrates MPAS dynamical core with CAM-SIMA by
+!> implementing the necessary APIs and managing their interaction.
+!>
+!> It implements bidirectional coupling between dynamics and physics states.
 module dyn_coupling
     ! Module(s) from CESM Share.
     use shr_kind_mod, only: kind_r8 => shr_kind_r8, &
@@ -663,7 +670,7 @@ pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t)
             ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
             ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
             !
-            ! In short, solve the below equation set for $T$ in terms of $\theta$, $\rho_d$ and $q_v$:
+            ! In all, solve the below equation set for $T$ in terms of $\theta$, $\rho_d$ and $q_v$:
             ! \begin{equation*}
             !     \begin{cases}
             !         \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
@@ -701,7 +708,7 @@ pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta)
             ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
             ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
             !
-            ! In short, solve the below equation set for $\theta$ in terms of $T$, $\rho_d$ and $q_v$:
+            ! In all, solve the below equation set for $\theta$ in terms of $T$, $\rho_d$ and $q_v$:
             ! \begin{equation*}
             !     \begin{cases}
             !         \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90
index 8a21b2a7..cfd9216c 100644
--- a/src/dynamics/mpas/dyn_grid.F90
+++ b/src/dynamics/mpas/dyn_grid.F90
@@ -1,3 +1,11 @@
+! SPDX-FileCopyrightText: Copyright (C) 2024 U.S. National Science Foundation National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+!> This module, part of the MPAS interface, integrates MPAS dynamical core with CAM-SIMA by
+!> implementing the necessary APIs and managing their interaction.
+!>
+!> It reads and uses information from MPAS mesh to initialize various model grids (e.g., dynamics,
+!> physics) for CAM-SIMA in terms of dynamics decomposition.
 module dyn_grid
     ! Module(s) from CAM-SIMA.
     use cam_grid_support, only: max_hcoordname_len
@@ -247,9 +255,9 @@ subroutine init_physics_grid()
         integer :: i
         integer :: ierr
         integer, pointer :: indextocellid(:)        ! Global indexes of cell centers.
-        real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (square meters).
-        real(kind_dyn_mpas), pointer :: latcell(:)  ! Cell center latitudes (radians).
-        real(kind_dyn_mpas), pointer :: loncell(:)  ! Cell center longitudes (radians).
+        real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (m2).
+        real(kind_dyn_mpas), pointer :: latcell(:)  ! Cell center latitudes (rad).
+        real(kind_dyn_mpas), pointer :: loncell(:)  ! Cell center longitudes (rad).
         type(physics_column_t), allocatable :: dyn_column(:) ! Grid and mapping information between global and local indexes.
 
         call dyn_debug_print(debugout_debug, subname // ' entered')
@@ -344,13 +352,13 @@ subroutine define_cam_grid()
         integer, pointer :: indextocellid(:)   ! Global indexes of cell centers.
         integer, pointer :: indextoedgeid(:)   ! Global indexes of edge nodes.
         integer, pointer :: indextovertexid(:) ! Global indexes of vertex nodes.
-        real(kind_dyn_mpas), pointer :: areacell(:)  ! Cell areas (square meters).
-        real(kind_dyn_mpas), pointer :: latcell(:)   ! Cell center latitudes (radians).
-        real(kind_dyn_mpas), pointer :: latedge(:)   ! Edge node latitudes (radians).
-        real(kind_dyn_mpas), pointer :: latvertex(:) ! Vertex node latitudes (radians).
-        real(kind_dyn_mpas), pointer :: loncell(:)   ! Cell center longitudes (radians).
-        real(kind_dyn_mpas), pointer :: lonedge(:)   ! Edge node longitudes (radians).
-        real(kind_dyn_mpas), pointer :: lonvertex(:) ! Vertex node longitudes (radians).
+        real(kind_dyn_mpas), pointer :: areacell(:)  ! Cell areas (m2).
+        real(kind_dyn_mpas), pointer :: latcell(:)   ! Cell center latitudes (rad).
+        real(kind_dyn_mpas), pointer :: latedge(:)   ! Edge node latitudes (rad).
+        real(kind_dyn_mpas), pointer :: latvertex(:) ! Vertex node latitudes (rad).
+        real(kind_dyn_mpas), pointer :: loncell(:)   ! Cell center longitudes (rad).
+        real(kind_dyn_mpas), pointer :: lonedge(:)   ! Edge node longitudes (rad).
+        real(kind_dyn_mpas), pointer :: lonvertex(:) ! Vertex node longitudes (rad).
 
         ! Global grid indexes. CAN be safely deallocated because its values are copied internally by
         ! `cam_grid_attribute_register` and `horiz_coord_create`.
@@ -360,7 +368,7 @@ subroutine define_cam_grid()
         ! just uses pointers to point at it internally.
         ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
         integer(kind_imap), pointer :: global_grid_map(:, :)
-        ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_register`
+        ! Cell areas (m2). CANNOT be safely deallocated because `cam_grid_attribute_register`
         ! just uses pointers to point at it internally.
         real(kind_r8), pointer :: cell_area(:)
         ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_register`
diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90
index dfaf2c18..62d1b6ce 100644
--- a/src/dynamics/mpas/stepon.F90
+++ b/src/dynamics/mpas/stepon.F90
@@ -1,3 +1,10 @@
+! SPDX-FileCopyrightText: Copyright (C) 2024 U.S. National Science Foundation National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+!> This module, part of the MPAS interface, integrates MPAS dynamical core with CAM-SIMA by
+!> implementing the necessary APIs and managing their interaction.
+!>
+!> It is a thin wrapper for plugging MPAS dynamical core into CAM-SIMA.
 module stepon
     ! Module(s) from CCPP.
     use ccpp_kinds, only: kind_phys