diff --git a/.dockerignore b/.dockerignore
index 267cdb79..12b85fe8 100644
--- a/.dockerignore
+++ b/.dockerignore
@@ -9,3 +9,6 @@
 !test/
 !.config_files.xml
 !docker
+!bin/
+!.lib/
+!.gitmodules
diff --git a/.github/scripts/branch_pr_issue_closer.py b/.github/scripts/branch_pr_issue_closer.py
index 1ad48ebe..71bb5235 100755
--- a/.github/scripts/branch_pr_issue_closer.py
+++ b/.github/scripts/branch_pr_issue_closer.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 
 """
 Script name:  branch_PR_issue_closer.py
@@ -21,52 +21,40 @@
 
 import re
 import sys
-import subprocess
-import shlex
 import argparse
 
 from github import Github
+from github import Auth
+
+###############
+#REGEX PATTERNS
+###############
+
+#Issue-closing Keywords are:
+#close, closes, closed
+#fix, fixes, fixed
+#resolve, resolves, resolved
+
+#The keywords are designed to match
+#the keywords that exist in Github
+#already for default branches, which
+#can be found here:
+#https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue
+
+#Create relevant regex patterns:
+_CLOSE_KEY              = r'close[sd]?'
+_FIX_KEY                = r'fix(e[sd])?'
+_RESOLVE_KEY            = r'resolve[sd]?'
+_KEYWORDS               = rf'{_CLOSE_KEY}|{_FIX_KEY}|{_RESOLVE_KEY}'
+_KEYWORDS_CAPTURE_GROUP = rf'(?P<keyword>{_KEYWORDS})'
+_ID_NUMBER              = r'\d+'
+_ID_CAPTURE_GROUP       = rf'(?P<id>{_ID_NUMBER})'
+_LINKED_ISSUE_PATTERN   = rf'{_KEYWORDS_CAPTURE_GROUP}\s*#{_ID_CAPTURE_GROUP}'
 
 #################
 #HELPER FUNCTIONS
 #################
 
-#+++++++++++++++++++++++++++++++++++++++++
-#Curl command needed to move project cards
-#+++++++++++++++++++++++++++++++++++++++++
-
-def  project_card_move(oa_token, column_id, card_id):
-
-    """
-    Currently pyGithub doesn't contain the methods required
-    to move project cards from one column to another, so
-    the unix curl command must be called directly, which is
-    what this function does.
-
-    The specific command-line call made is:
-
-    curl -H "Authorization: token OA_token" -H \
-    "Accept: application/vnd.github.inertia-preview+json" \
-    -X POST -d '{"position":"top", "column_id":<column_id>}' \
-    https://api.github.com/projects/columns/cards/<card_id>/moves
-
-    """
-
-    #create required argument strings from inputs:
-    github_oa_header = f''' "Authorization: token {oa_token}" '''
-    github_url_str = f'''https://api.github.com/projects/columns/cards/{card_id}/moves'''
-    json_post_inputs = f''' '{{"position":"top", "column_id":{column_id}}}' '''
-
-    #Create curl command line string:
-    curl_cmdline = '''curl -H '''+github_oa_header+''' -H "Accept: application/vnd.github.inertia-preview+json" -X POST -d '''+\
-                   json_post_inputs+''' '''+github_url_str
-
-    #Split command line string into argument list:
-    curl_arg_list = shlex.split(curl_cmdline)
-
-    #Run command using subprocess:
-    subprocess.run(curl_arg_list, check=True)
-
 #++++++++++++++++++++++++++++++
 #Input Argument parser function
 #++++++++++++++++++++++++++++++
@@ -135,176 +123,129 @@ def _main_prog():
     #Log-in to github API using token
     #++++++++++++++++++++++++++++++++
 
-    ghub = Github(token)
+    auth = Auth.Token(token)
+    ghub = Github(auth=auth)
 
-    #+++++++++++++++++++++
+    #+++++++++++++++++++++++++
     #Open ESCOMP/CAM-SIMA repo
-    #+++++++++++++++++++++
+    #+++++++++++++++++++++++++
 
     cam_repo = ghub.get_repo("ESCOMP/CAM-SIMA")
 
-    #+++++++++++++++++++++++++++++
-    #Get triggering commit message
-    #+++++++++++++++++++++++++++++
+    #++++++++++++++++++++++++++++++
+    #Get PRs associated with commit
+    #++++++++++++++++++++++++++++++
 
     github_commit = cam_repo.get_commit(trigger_sha)
 
-    commit_message = github_commit.commit.message
+    commit_prs = github_commit.get_pulls()
 
-    #+++++++++++++++++++++++++++++++
-    #Search for github PR merge text
-    #+++++++++++++++++++++++++++++++
+    pr_nums = [pr.number for pr in commit_prs]
 
-    #Compile Pull Request merge text expression:
-    pr_merge_pattern = re.compile(r'Merge pull request ')
 
-    #Search for merge text, starting at beginning of message:
-    commit_msg_match = pr_merge_pattern.match(commit_message)
+    #If list is empty, then no PRs are associated
+    #with this commit, so go ahead and close:
+    if not pr_nums:
+        endmsg = f"No PRs associated with commit:\n{trigger_sha}\n"
+        endmsg += " so issue-closing script is stopping here."
+        end_script(endmsg)
 
-    #Check if match exists:
-    if commit_msg_match is not None:
-        #If it does then pull out text immediately after message:
-        post_msg_text = commit_message[commit_msg_match.end():]
+    #++++++++++++++++++++++++++++
+    #Loop over all associated PRs
+    #++++++++++++++++++++++++++++
 
-        #Split text into individual words:
-        post_msg_word_list = post_msg_text.split()
+    for pr_num in pr_nums:
 
-        #Extract first word:
-        first_word = post_msg_word_list[0]
+        #+++++++++++++++++++++++++++++++++++++
+        #Check that PR has in fact been merged
+        #+++++++++++++++++++++++++++++++++++++
 
-        #Print merged pr number to screen:
-        print(f"Merged PR: {first_word}")
+        #Extract pull request info:
+        merged_pull = cam_repo.get_pull(pr_num)
 
-        try:
-            #Try assuming the word is just a number:
-            pr_num = int(first_word[1:]) #ignore "#" symbol
-        except ValueError:
-            #If the conversion fails, then this is likely not a real PR merge, so end the script:
-            endmsg = "No Pull Request number was found in the commit message, so there is nothing for the script to do."
+        #If pull request has not been merged, then exit script:
+        if not merged_pull.merged:
+            endmsg = f"Pull request #{pr_num} associated with commit:\n{trigger_sha}\n"
+            endmsg += "was not actually merged, so the script will not close anything."
             end_script(endmsg)
 
-    else:
-        endmsg = "This push commit does not appear to be a merged pull request, so the script will do nothing."
-        end_script(endmsg)
+        #++++++++++++++++++++++++++++++++++++++++
+        #Check that PR was not for default branch
+        #++++++++++++++++++++++++++++++++++++++++
 
-    #+++++++++++++++++++++++++++++++++++++
-    #Check that PR has in fact been merged
-    #+++++++++++++++++++++++++++++++++++++
+        #Determine default branch on repo:
+        default_branch = cam_repo.default_branch
 
-    #Extract pull request info:
-    merged_pull = cam_repo.get_pull(pr_num)
+        #Extract merged branch from latest Pull request:
+        merged_branch = merged_pull.base.ref
 
-    #If pull request has not been merged, then exit script:
-    if not merged_pull.merged:
-        endmsg = "Pull request in commit message was not actually merged, so the script will not close anything."
-        end_script(endmsg)
+        #If PR was to default branch, then exit script (as github will handle it automatically):
+        if merged_branch == default_branch:
+            endmsg = f"Pull request #{pr_num} was merged into default repo branch. "
+            endmsg += "Thus issue is closed automatically"
+            end_script(endmsg)
 
-    #++++++++++++++++++++++++++++++++++++++++
-    #Check that PR was not for default branch
-    #++++++++++++++++++++++++++++++++++++++++
+        #++++++++++++++++++++++++++++++++++++++
+        #Create integer list of all open issues:
+        #++++++++++++++++++++++++++++++++++++++
 
-    #Determine default branch on repo:
-    default_branch = cam_repo.default_branch
+        #Extract list of open issues from repo:
+        open_repo_issues = cam_repo.get_issues(state='open')
 
-    #Extract merged branch from latest Pull request:
-    merged_branch = merged_pull.base.ref
+        #Collect all open repo issues:
+        open_issues = [issue.number for issue in open_repo_issues]
 
-    #If PR was to default branch, then exit script (as github will handle it automatically):
-    if merged_branch == default_branch:
-        endmsg = "Pull request ws merged into default repo branch. Thus issue is closed automatically"
-        end_script(endmsg)
+        #+++++++++++++++++++++++++++++++++++++++++++++
+        #Create integer list of all open pull requests
+        #+++++++++++++++++++++++++++++++++++++++++++++
 
-    #++++++++++++++++++++++++++++++++++++++
-    #Create integer list of all open issues:
-    #++++++++++++++++++++++++++++++++++++++
+        #Extract list of open PRs from repo:
+        open_repo_pulls = cam_repo.get_pulls(state='open')
 
-    #Extract list of open issues from repo:
-    open_repo_issues = cam_repo.get_issues(state='open')
+        #Collect all open pull requests:
+        open_pulls = [pr.number for pr in open_repo_pulls]
 
-    #Collect all open repo issues:
-    open_issues = [issue.number for issue in open_repo_issues]
+        #+++++++++++++++++++++++++++++++++++++++++++++++++
+        #Check if one of the keywords exists in PR message
+        #+++++++++++++++++++++++++++++++++++++++++++++++++
 
-    #+++++++++++++++++++++++++++++++++++++++++++++
-    #Create integer list of all open pull requests
-    #+++++++++++++++++++++++++++++++++++++++++++++
+        #Compile regex patterns into object:
+        keyword_pattern = re.compile(_LINKED_ISSUE_PATTERN)
 
-    #Extract list of open PRs from repo:
-    open_repo_pulls = cam_repo.get_pulls(state='open')
+        #Extract (lower case) Pull Request message:
+        pr_msg_lower = merged_pull.body.lower()
 
-    #Collect all open pull requests:
-    open_pulls = [pr.number for pr in open_repo_pulls]
+        #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        #Extract issue and PR numbers associated with found keywords in merged PR message
+        #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-    #+++++++++++++++++++++++++++++++++++++++++++++++++
-    #Check if one of the keywords exists in PR message
-    #+++++++++++++++++++++++++++++++++++++++++++++++++
+        #Create new "close" issues list:
+        close_issues = []
 
-    #Keywords are:
-    #close, closes, closed
-    #fix, fixes, fixed
-    #resolve, resolves, resolved
+        #Create new "closed" PR list:
+        close_pulls = []
 
-    #Create regex pattern to find keywords:
-    keyword_pattern = re.compile(r'(^|\s)close(\s|s\s|d\s)|(^|\s)fix(\s|es\s|ed\s)|(^|\s)resolve(\s|s\s|d\s)')
+        #Create iterator of all keyword/id pairs:
+        word_matches = keyword_pattern.finditer(pr_msg_lower, re.IGNORECASE)
 
-    #Extract (lower case) Pull Request message:
-    pr_msg_lower = merged_pull.body.lower()
+        #Go through all matches to pull out PR and issue numbers:
+        found_ids = set()
+        for match in word_matches:
+            issue_dict = match.groupdict()
+            issue_num  = int(issue_dict['id'].lstrip('0'))
+            found_ids.add(issue_num)
 
-    #search for at least one keyword:
-    if keyword_pattern.search(pr_msg_lower) is not None:
-        #If at least one keyword is found, then determine location of every keyword instance:
-        word_matches = keyword_pattern.finditer(pr_msg_lower)
-    else:
-        endmsg = "Pull request was merged without using any of the keywords.  Thus there are no issues to close."
-        end_script(endmsg)
+        #End script if no keyword/id pairs were found:
+        if not found_ids:
+            endmsg = f"Pull request #{pr_num} was merged without using any of the keywords. "
+            endmsg += "Thus there are no issues to close."
+            end_script(endmsg)
+
+        close_pulls = list(found_ids.intersection(open_pulls))
+        close_issues = list(found_ids.intersection(open_issues))
 
-    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    #Extract issue and PR numbers associated with found keywords in merged PR message
-    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-    #create issue pattern ("the number symbol {#} + a number"),
-    #which ends with either a space, a comma, a period, or
-    #the end of the string itself:
-    issue_pattern = re.compile(r'#[0-9]+(\s|,|$)|.')
-
-    #Create new "close" issues list:
-    close_issues = []
-
-    #Create new "closed" PR list:
-    close_pulls = []
-
-    #Search text right after keywords for possible issue numbers:
-    for match in word_matches:
-
-        #create temporary string starting at end of match:
-        tmp_msg_str = pr_msg_lower[match.end():]
-
-        #Check if first word matches issue pattern:
-        if issue_pattern.match(tmp_msg_str) is not None:
-
-            #If so, then look for an issue number immediately following
-            first_word = tmp_msg_str.split()[0]
-
-            #Extract issue number from first word:
-            try:
-                #First try assuming the string is just a number
-                issue_num = int(first_word[1:]) #ignore "#" symbol
-            except ValueError:
-                #If not, then ignore last letter:
-                try:
-                    issue_num = int(first_word[1:-1])
-                except ValueError:
-                    #If ignoring the first and last letter doesn't work,
-                    #then the match was likely a false positive,
-                    #so set the issue number to one that will never be found:
-                    issue_num = -9999
-
-            #Check if number is actually for a PR (as opposed to an issue):
-            if issue_num in open_pulls:
-                #Add PR number to "close pulls" list:
-                close_pulls.append(issue_num)
-            elif issue_num in open_issues:
-                #If in fact an issue, then add to "close issues" list:
-                close_issues.append(issue_num)
+
+    #+++END REFERENCED PR LOOP+++
 
     #If no issue numbers are present after any of the keywords, then exit script:
     if not close_issues and not close_pulls:
@@ -321,177 +262,19 @@ def _main_prog():
         print("PRs referenced by the merged PR: "+", ".join(\
               str(pull) for pull in close_pulls))
 
-    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    #Determine name of project associated with merged Pull Request
-    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-    #Pull-out all projects from repo:
-    projects = cam_repo.get_projects()
-
-    #Initalize modified project name:
-    proj_mod_name = None
-
-    #Loop over all repo projects:
-    for project in projects:
-
-        #Pull-out columns from each project:
-        proj_columns = project.get_columns()
-
-        #Loop over columns:
-        for column in proj_columns:
-
-            #check if column name is "Completed tags"
-            if column.name == "Completed tags":
-                #If so, then extract cards:
-                cards = column.get_cards()
-
-                #Loop over cards:
-                for card in cards:
-                    #Extract card content:
-                    card_content = card.get_content()
-
-                    #Next, check if card number exists and matches merged PR number:
-                    if card_content is not None and card_content.number == pr_num:
-                        #If so, and if Project name is None, then set string:
-                        if proj_mod_name is None:
-                            proj_mod_name = project.name
-                            #Break out of card loop:
-                            break
-
-                        #If already set, then somehow merged PR is in two different projects,
-                        #which is not what this script is expecting, so just exit:
-                        endmsg = "Merged Pull Request found in two different projects, so script will do nothing."
-                        end_script(endmsg)
-
-    #Print project name associated with merged PR:
-    print(f"merged PR project name: {proj_mod_name}")
-
-    #++++++++++++++++++++++++++++++++++++++++
-    #Extract repo project "To do" card issues
-    #++++++++++++++++++++++++++++++++++++++++
-
-    #Initalize issue counting dictionary:
-    proj_issues_count = {}
-
-    #Initalize issue id to project card id dictionary:
-    proj_issue_card_ids = {}
-
-    #Initialize list for issues that have already been closed:
-    already_closed_issues = []
-
-    #Loop over all repo projects:
-    for project in projects:
-
-        #Next, pull-out columns from each project:
-        proj_columns = project.get_columns()
-
-        #Loop over columns:
-        for column in proj_columns:
-            #Check if column name is "To do"
-            if column.name == "To do":
-                #If so, then extract cards:
-                cards = column.get_cards()
-
-                #Loop over cards:
-                for card in cards:
-                    #Extract card content:
-                    card_content = card.get_content()
-
-                    #Next, check if card issue number matches any of the "close" issue numbers from the PR:
-                    if card_content is not None and card_content.number in close_issues:
-
-                        #If so, then check if issue number is already in proj_issues_count:
-                        if card_content.number in proj_issues_count:
-                            #Add one to project issue counter:
-                            proj_issues_count[card_content.number] += 1
-
-                            #Also add issue id and card id to id dictionary used for card move, if in relevant project:
-                            if project.name == proj_mod_name:
-                                proj_issue_card_ids[card_content.number] = card.id
-
-                        else:
-                            #If not, then append to project issues count dictionary:
-                            proj_issues_count[card_content.number] = 1
-
-                            #Also add issue id and card id to id dictionary used for card move, if in relevant project:
-                            if project.name == proj_mod_name:
-                                proj_issue_card_ids[card_content.number] = card.id
-
-            #Otherwise, check if column name matches "closed issues" column:
-            elif column.name == "closed issues" and project.name == proj_mod_name:
-                #Save column id:
-                column_target_id = column.id
-
-                #Extract cards:
-                closed_cards = column.get_cards()
-
-                #Loop over cards:
-                for closed_card in closed_cards:
-                    #Extract card content:
-                    closed_card_content = closed_card.get_content()
-
-                    #Check if card issue number matches any of the "close" issue numbers from the PR:
-                    if closed_card_content is not None and closed_card_content.number in close_issues:
-                        #If issue number matches, then it likely means the same
-                        #commit message or issue number reference was used in multiple
-                        #pushes to the same repo (e.g., for a PR and then a tag). Thus
-                        #the issue should be marked as "already closed":
-                        already_closed_issues.append(closed_card_content.number)
-
-    #Remove all issues from issue dictionary that are "already closed":
-    for already_closed_issue_num in already_closed_issues:
-        if already_closed_issue_num in proj_issues_count:
-            proj_issues_count.pop(already_closed_issue_num)
-
-    #If no project cards are found that match the issue, then exit script:
-    if not proj_issues_count:
-        endmsg = "No project cards match the issue being closed, so the script will do nothing."
-        end_script(endmsg)
+    #++++++++++++++++++++++++++++++++++++++++++++++
+    #Attempt to close all referenced issues and PRs
+    #++++++++++++++++++++++++++++++++++++++++++++++
 
-    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    #Check if the number of "To-do" project cards matches the total number
-    #of merged PRs for each 'close' issue.
-    #
-    #Then, close all issues for which project cards equals merged PRs
-    #
-    #If not, then simply move the project card to the relevant project's
-    #"closed issues" column.
-    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-    #Loop over project issues and counts that have been "closed" by merged PR:
-    for issue_num, issue_count in proj_issues_count.items():
-
-        #If issue count is just one, then close issue:
-        if issue_count == 1:
-            #Extract github issue object:
-            cam_issue = cam_repo.get_issue(number=issue_num)
-            #Close issue:
-            cam_issue.edit(state='closed')
-            print(f"Issue #{issue_num} has been closed.")
-        else:
-            #Extract card id from id dictionary:
-            if issue_num in proj_issue_card_ids:
-                card_id = proj_issue_card_ids[issue_num]
-            else:
-                #If issue isn't in dictionary, then it means the issue
-                #number was never found in the "To do" column, which
-                #likely means the user either referenced the wrong
-                #issue number, or the issue was never assigned to the
-                #project.  Warn user and then exit with a non-zero
-                #error so that the Action fails:
-                endmsg = 'Issue #{} was not found in the "To Do" Column of the "{}" project.\n' \
-                         'Either the wrong issue number was referenced, or the issue was never ' \
-                         'attached to the project.'.format(issue_num, proj_mod_name)
-                print(endmsg)
-                sys.exit(1)
-
-            #Then move the card on the relevant project page to the "closed issues" column:
-            project_card_move(token.strip(), column_target_id, card_id)
-
-    #++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    #Finally, close all Pull Requests in "close_pulls" list:
-    #++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    #Loop over referenced issues:
+    for issue_num in close_issues:
+        #Extract github issue object:
+        cam_issue = cam_repo.get_issue(number=issue_num)
+        #Close issue:
+        cam_issue.edit(state='closed')
+        print(f"Issue #{issue_num} has been closed.")
 
+    #Loop over referenced PRs:
     for pull_num in close_pulls:
         #Extract Pull request object:
         cam_pull = cam_repo.get_pull(number=pull_num)
diff --git a/.github/workflows/branch_push_workflow.yml b/.github/workflows/branch_push_workflow.yml
index 94f4414a..d06af346 100644
--- a/.github/workflows/branch_push_workflow.yml
+++ b/.github/workflows/branch_push_workflow.yml
@@ -20,12 +20,12 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     # Acquire github action routines
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     # Acquire specific version of python
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v4
+    - name: Set up Python 3.11
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.10' # Semantic version range syntax or exact version of a Python version
+        python-version: '3.11' # Semantic version range syntax or exact version of a Python version
     # Install required python packages
     - name: Install dependencies
       run: |
diff --git a/.github/workflows/fleximod_test.yaml b/.github/workflows/fleximod_test.yaml
new file mode 100644
index 00000000..36e68b19
--- /dev/null
+++ b/.github/workflows/fleximod_test.yaml
@@ -0,0 +1,18 @@
+on:
+  pull_request:
+jobs:
+  fleximod-test:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        # oldest supported and latest supported
+        python-version: ["3.8", "3.x"]
+    steps:
+      - id: checkout-CESM
+        uses: actions/checkout@v4
+      - id: run-fleximod
+        run: |
+          $GITHUB_WORKSPACE/bin/git-fleximod update
+          $GITHUB_WORKSPACE/bin/git-fleximod test
+        
+          
diff --git a/.gitmodules b/.gitmodules
index c0b72575..2c426e6e 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,7 +1,7 @@
 [submodule "ccpp-framework"]
         path          = ccpp_framework
         url           = https://github.com/NCAR/ccpp-framework
-        fxtag         = 2024-07-11-dev
+        fxtag         = 2024-07-19-dev
         fxrequired    = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/NCAR/ccpp-framework
 [submodule "history"]
@@ -20,7 +20,7 @@
 [submodule "ncar-physics"]
         path          = src/physics/ncar_ccpp
 	url           = https://github.com/peverwhee/atmospheric_physics
-	fxtag         = diagnostics
+	fxtag         = 61bd9d3dac2abb11bb1e44a2ca34b401da0f44b1
         fxrequired    = AlwaysRequired
         fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
 [submodule "ccs_config"]
diff --git a/.lib/git-fleximod/git_fleximod/cli.py b/.lib/git-fleximod/git_fleximod/cli.py
index fdcf102a..b6f728f8 100644
--- a/.lib/git-fleximod/git_fleximod/cli.py
+++ b/.lib/git-fleximod/git_fleximod/cli.py
@@ -2,7 +2,7 @@
 import argparse
 from git_fleximod import utils
 
-__version__ = "0.8.2"
+__version__ = "0.8.4"
 
 def find_root_dir(filename=".gitmodules"):
     """ finds the highest directory in tree
diff --git a/.lib/git-fleximod/git_fleximod/git_fleximod.py b/.lib/git-fleximod/git_fleximod/git_fleximod.py
index 4595cd2a..50e0ef83 100755
--- a/.lib/git-fleximod/git_fleximod/git_fleximod.py
+++ b/.lib/git-fleximod/git_fleximod/git_fleximod.py
@@ -195,18 +195,19 @@ def submodules_status(gitmodules, root_dir, toplevel=False, depth=0):
         submod = init_submodule_from_gitmodules(gitmodules, name, root_dir, logger)
             
         result,n,l,t = submod.status()
-        testfails += t
-        localmods += l
-        needsupdate += n
         if toplevel or not submod.toplevel():
             print(wrapper.fill(result))
-        subdir = os.path.join(root_dir, submod.path)
-        if os.path.exists(os.path.join(subdir, ".gitmodules")):
-            submod = GitModules(logger, confpath=subdir)
-            t,l,n = submodules_status(submod, subdir, depth=depth+1)
             testfails += t
             localmods += l
             needsupdate += n
+        subdir = os.path.join(root_dir, submod.path)
+        if os.path.exists(os.path.join(subdir, ".gitmodules")):
+            gsubmod = GitModules(logger, confpath=subdir)
+            t,l,n = submodules_status(gsubmod, subdir, depth=depth+1)
+            if toplevel or not submod.toplevel():
+                testfails += t
+                localmods += l
+                needsupdate += n
             
     return testfails, localmods, needsupdate
 
diff --git a/.lib/git-fleximod/git_fleximod/submodule.py b/.lib/git-fleximod/git_fleximod/submodule.py
index 48657272..70a3018a 100644
--- a/.lib/git-fleximod/git_fleximod/submodule.py
+++ b/.lib/git-fleximod/git_fleximod/submodule.py
@@ -92,7 +92,7 @@ def status(self):
                 needsupdate = True
             else:
                 result = f"e {self.name:>20} has no fxtag defined in .gitmodules{optional}"
-                testfails = True
+                testfails = False
         else:
             with utils.pushd(smpath):
                 git = GitInterface(smpath, self.logger)
@@ -103,15 +103,23 @@ def status(self):
                     needsupdate = True
                     return result, needsupdate, localmods, testfails                    
                 rurl = git.git_operation("ls-remote","--get-url").rstrip()
-                line = git.git_operation("log", "--pretty=format:\"%h %d").partition('\n')[0]
+                line = git.git_operation("log", "--pretty=format:\"%h %d\"").partition('\n')[0]
                 parts = line.split()
                 ahash = parts[0][1:]
+                atag = None
                 if len(parts) > 3:
-                    atag = parts[3][:-1]
-                else:
-                    atag = None
+                    idx = 0
+                    while idx < len(parts)-1:
+                        idx = idx+1
+                        if parts[idx] == 'tag:':
+                            atag = parts[idx+1]
+                            while atag.endswith(')') or atag.endswith(',') or atag.endswith("\""):
+                                atag = atag[:-1]
+                            if atag == self.fxtag:
+                                break
+
                 
-                #print(f"line is {line} ahash is {ahash} atag is {atag}")
+                #print(f"line is {line} ahash is {ahash} atag is {atag} {parts}")
                 #                atag = git.git_operation("describe", "--tags", "--always").rstrip()
                 # ahash =  git.git_operation("rev-list", "HEAD").partition("\n")[0]
                     
@@ -122,9 +130,11 @@ def status(self):
                 if self.fxtag and atag == self.fxtag:
                     result = f"  {self.name:>20} at tag {self.fxtag}"
                     recurse = True
-                elif self.fxtag and ahash[: len(self.fxtag)] == self.fxtag:
+                    testfails = False
+                elif self.fxtag and (ahash[: len(self.fxtag)] == self.fxtag or (self.fxtag.find(ahash)==0)):
                     result = f"  {self.name:>20} at hash {ahash}"
                     recurse = True
+                    testfails = False
                 elif atag == ahash:
                     result = f"  {self.name:>20} at hash {ahash}"
                     recurse = True
@@ -133,14 +143,17 @@ def status(self):
                     testfails = True
                     needsupdate = True
                 else:
-                    result = f"e {self.name:>20} has no fxtag defined in .gitmodules, module at {atag}"
-                    testfails = True
+                    if atag:
+                        result = f"e {self.name:>20} has no fxtag defined in .gitmodules, module at {atag}"
+                    else:
+                        result = f"e {self.name:>20} has no fxtag defined in .gitmodules, module at {ahash}"
+                    testfails = False
                     
                 status = git.git_operation("status", "--ignore-submodules", "-uno")
                 if "nothing to commit" not in status:
                     localmods = True
                     result = "M" + textwrap.indent(status, "                      ")
-        
+#        print(f"result {result} needsupdate {needsupdate} localmods {localmods} testfails {testfails}")
         return result, needsupdate, localmods, testfails
 
     
diff --git a/.lib/git-fleximod/pyproject.toml b/.lib/git-fleximod/pyproject.toml
index 9cff1423..850e57d5 100644
--- a/.lib/git-fleximod/pyproject.toml
+++ b/.lib/git-fleximod/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "git-fleximod"
-version = "0.8.2"
+version = "0.8.4"
 description = "Extended support for git-submodule and git-sparse-checkout"
 authors = ["Jim Edwards <jedwards@ucar.edu>"]
 maintainers = ["Jim Edwards <jedwards@ucar.edu>"]
diff --git a/.lib/git-fleximod/tbump.toml b/.lib/git-fleximod/tbump.toml
index b4eed7d4..bd82c557 100644
--- a/.lib/git-fleximod/tbump.toml
+++ b/.lib/git-fleximod/tbump.toml
@@ -2,7 +2,7 @@
 github_url = "https://github.com/jedwards4b/git-fleximod/"
 
 [version]
-current = "0.8.2"
+current = "0.8.4"
 
 # Example of a semver regexp.
 # Make sure this matches current_version before
diff --git a/ccpp_framework b/ccpp_framework
index 0f823272..5f338ddf 160000
--- a/ccpp_framework
+++ b/ccpp_framework
@@ -1 +1 @@
-Subproject commit 0f8232724975c13289cad390c9a71fa2c6a9bff4
+Subproject commit 5f338ddf02638c06548e54e0a218d154b34faff9
diff --git a/cime_config/cam_autogen.py b/cime_config/cam_autogen.py
index 0520d602..eeb31229 100644
--- a/cime_config/cam_autogen.py
+++ b/cime_config/cam_autogen.py
@@ -615,6 +615,11 @@ def generate_physics_suites(build_cache, preproc_defs, host_name,
         ufiles_str = datatable_report(cap_output_file, request, ";")
         utility_files = ufiles_str.split(';')
         _update_genccpp_dir(utility_files, genccpp_dir)
+        request = DatatableReport("dependencies")
+        dep_str = datatable_report(cap_output_file, request, ";")
+        if len(dep_str) > 0:
+            dependency_files = dep_str.split(';')
+            _update_genccpp_dir(dependency_files, genccpp_dir)
         ##XXgoldyXX: ^ Temporary fix:
     # End if
 
diff --git a/docker/Dockerfile b/docker/Dockerfile
index 1d6acc6c..098db37d 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -1,6 +1,6 @@
 # parts of CAM require x86 architecture (gptl, which relies on the rdtsc x86 assembly instruction)
 # esmf is am image you are expected to have built. Read the README file for instructions
-FROM esmf:latest
+FROM --platform=linux/amd64 esmf:latest
 
 ###################################################
 ## Install necessary packages
@@ -11,13 +11,18 @@ RUN dnf -y update \
       git \
       hostname \
       m4 \
-      python \
+      python39 \
+      pip \
       sudo \
       svn \
       tree \
       vim \
     && dnf clean all
 
+RUN ln -s $(which python3) /usr/bin/python && \
+    pip install --upgrade pip && \
+    pip install --upgrade setuptools
+
 ###################################################
 ## Make sure the mpi compilers can be found
 ###################################################
@@ -50,7 +55,7 @@ USER cam_sima_user
 WORKDIR /home/cam_sima_user/CAM-SIMA
 
 # pull the dependencies
-RUN ./manage_externals/checkout_externals
+RUN ./bin/git-fleximod update
 
 # Copy in the machine information for the container
 RUN cp /home/cam_sima_user/CAM-SIMA/docker/config_machines.xml /home/cam_sima_user/CAM-SIMA/ccs_config/machines/
@@ -79,7 +84,7 @@ RUN ./xmlchange STOP_N=5
 RUN chmod +x /home/cam_sima_user/CAM-SIMA/docker/ftp_download.sh
 RUN /home/cam_sima_user/CAM-SIMA/docker/ftp_download.sh
 
-# # add the snapshot file
+# add the snapshot file
 RUN echo "ncdata='/home/cam_sima_user/run_heldsuarez_cam6_nt2_bigg_try005.cam.h5.0001-01-01-00000.nc'" >> user_nl_cam
 
 RUN ./case.build
\ No newline at end of file
diff --git a/docker/Dockerfile.musica b/docker/Dockerfile.musica
index b02cf7b6..0f59f21d 100644
--- a/docker/Dockerfile.musica
+++ b/docker/Dockerfile.musica
@@ -1,6 +1,6 @@
 # parts of CAM require x86 architecture (gptl, which relies on the rdtsc x86 assembly instruction)
 # esmf is am image you are expected to have built. Read the README file for instructions
-FROM esmf:latest
+FROM --platform=linux/amd64 esmf:latest
 
 ###################################################
 ## Install necessary packages
@@ -11,13 +11,18 @@ RUN dnf -y update \
       git \
       hostname \
       m4 \
-      python \
+      python39 \
+      pip \
       sudo \
       svn \
       tree \
       vim \
     && dnf clean all
 
+RUN ln -s $(which python3) /usr/bin/python && \
+    pip install --upgrade pip && \
+    pip install --upgrade setuptools
+
 ###################################################
 ## Make sure the mpi compilers can be found
 ###################################################
@@ -36,37 +41,22 @@ RUN cd pnetcdf-1.12.3 && \
 
 ENV FC=gfortran
 
-###################################################
-## Build and install json-fortran
-###################################################
-RUN curl -LO https://github.com/jacobwilliams/json-fortran/archive/8.2.0.tar.gz \
-    && tar -zxvf 8.2.0.tar.gz \
-    && cd json-fortran-8.2.0 \
-    && mkdir build \
-    && cd build \
-    && cmake -D SKIP_DOC_GEN:BOOL=TRUE .. \
-    && make install -j 8
-
-# add a symlink
-RUN ln -s /usr/local/jsonfortran-gnu-8.2.0/lib/libjsonfortran.a /usr/local/lib/libjsonfortran.a
-
 ###################################################
 ## Build and install MUSICA
 ###################################################
 
-RUN git clone https://github.com/NCAR/musica.git
+RUN git clone https://github.com/NCAR/musica.git \
+    && cd musica \
+    && git checkout 2a5eeaac982a3eb80b96d1e2087b91b301d1e748
+
 RUN mkdir /musica/build \
     && cd /musica/build \
-    && export JSON_FORTRAN_HOME="/usr/local/jsonfortran-gnu-8.2.0" \
     && cmake \
              -D ENABLE_TESTS=OFF \
-             -D ENABLE_TUVX=OFF \
-          .. \
+             -D MUSICA_BUILD_FORTRAN_INTERFACE=ON \
+             .. \
     && make install -j 8
 
-# add a symlink
-RUN ln -s /usr/local/musica-0.3.0/lib64/libmusica.a /usr/local/lib/libmusica.a
-
 ###################################################
 ## Build CAM-SIMA
 ###################################################
@@ -83,7 +73,7 @@ USER cam_sima_user
 WORKDIR /home/cam_sima_user/CAM-SIMA
 
 # pull the dependencies
-RUN ./manage_externals/checkout_externals
+RUN ./bin/git-fleximod update
 
 # Copy in the machine information for the container
 RUN cp /home/cam_sima_user/CAM-SIMA/docker/config_machines.xml /home/cam_sima_user/CAM-SIMA/ccs_config/machines/
@@ -104,7 +94,7 @@ WORKDIR $CASE_NAME
 RUN ./case.setup
 
 RUN ./xmlchange CAM_CONFIG_OPTS="--dyn none --physics-suites musica"
-RUN ./xmlchange CAM_LINKED_LIBS="-lmusica -ljsonfortran"
+RUN ./xmlchange CAM_LINKED_LIBS="-lmusica-fortran -lmusica -lyaml-cpp"
 RUN ./xmlchange ROF_NCPL=48
 RUN ./xmlchange STOP_OPTION=nsteps
 RUN ./xmlchange STOP_N=5
diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90
index cedbea8b..fff0f506 100644
--- a/src/control/cam_comp.F90
+++ b/src/control/cam_comp.F90
@@ -546,7 +546,6 @@ subroutine cam_register_constituents(cam_runtime_opts)
       integer                                        :: errflg
       character(len=512)                             :: errmsg
       type(ccpp_constituent_prop_ptr_t), pointer     :: const_props(:)
-      type(ccpp_constituent_properties_t), allocatable, target :: dynamic_constituents(:)
       character(len=*), parameter :: subname = 'cam_register_constituents: '
 
       ! Initalize error flag and message:
diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/Makefile.in.CESM
index a0bff7ed..f4016f0e 100644
--- a/src/dynamics/mpas/Makefile.in.CESM
+++ b/src/dynamics/mpas/Makefile.in.CESM
@@ -130,5 +130,5 @@ externals: $(AUTOCLEAN_DEPS)
 subdrv: driver/dyn_mpas_subdriver.o
 
 %.o: %.F90 dycore frame ops
-	( cd $(<D); $(FC) $(CPPFLAGS) $(FFLAGS) -c $(<F) $(CPPINCLUDES) $(FCINCLUDES) -I../framework -I../operators -I../core_$(CORE) )
+	( cd $(<D); $(FC) $(CPPFLAGS) $(FFLAGS) -c $(<F) $(CPPINCLUDES) $(FCINCLUDES) -I../framework -I../operators -I../core_$(CORE) -I../core_$(CORE)/dynamics )
 	$(LN) $(@) .
diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
index 9d2085f1..b21cdcfe 100644
--- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
+++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
@@ -11,16 +11,25 @@ module dyn_mpas_subdriver
 
     ! Modules from external libraries.
     use mpi, only: mpi_comm_null, mpi_comm_rank, mpi_success
-    use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, &
-                   pio_inq_varid, pio_inq_varndims, pio_noerr
+    use pio, only: pio_char, pio_int, pio_real, pio_double, &
+                   file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, &
+                   pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr
 
     ! Modules from MPAS.
+    use atm_core, only: atm_mpas_init_block
     use atm_core_interface, only: atm_setup_core, atm_setup_domain
+    use atm_time_integration, only: mpas_atm_dynamics_init
+    use mpas_atm_dimensions, only: mpas_atm_set_dims
+    use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group
+    use mpas_atm_threading, only: mpas_atm_threading_init
+    use mpas_attlist, only: mpas_modify_att
     use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2
+    use mpas_constants, only: mpas_constants_compute_derived
     use mpas_derived_types, only: core_type, domain_type, &
                                   mpas_pool_type, mpas_pool_field_info_type, &
                                   mpas_pool_character, mpas_pool_real, mpas_pool_integer, &
                                   mpas_stream_type, mpas_stream_noerr, &
+                                  mpas_time_type, mpas_start_time, &
                                   mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, &
                                   field0dchar, field1dchar, &
                                   field0dinteger, field1dinteger, field2dinteger, field3dinteger, &
@@ -31,13 +40,17 @@ module dyn_mpas_subdriver
     use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2
     use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, &
                                mpas_readstream, mpas_writestream, mpas_writestreamatt
-    use mpas_kind_types, only: rkind, strkind
+    use mpas_kind_types, only: rkind, r4kind, r8kind, strkind
     use mpas_pool_routines, only: mpas_pool_create_pool, mpas_pool_destroy_pool, mpas_pool_get_subpool, &
                                   mpas_pool_add_config, mpas_pool_get_config, &
                                   mpas_pool_get_array, &
-                                  mpas_pool_get_dimension, &
-                                  mpas_pool_get_field, mpas_pool_get_field_info
+                                  mpas_pool_add_dimension, mpas_pool_get_dimension, &
+                                  mpas_pool_get_field, mpas_pool_get_field_info, &
+                                  mpas_pool_initialize_time_levels
     use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex
+    use mpas_string_utils, only: mpas_string_replace
+    use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time
+    use mpas_vector_operations, only: mpas_initialize_vectors
 
     implicit none
 
@@ -61,6 +74,7 @@ end subroutine model_error_if
 
         logical, public :: debug_output = .false.
 
+        ! Initialized by `dyn_mpas_init_phase1`.
         integer :: log_unit = output_unit
         integer :: mpi_comm = mpi_comm_null
         integer :: mpi_rank = 0
@@ -71,6 +85,18 @@ end subroutine model_error_if
 
         type(core_type), pointer :: corelist => null()
         type(domain_type), pointer :: domain_ptr => null()
+
+        ! Initialized by `dyn_mpas_init_phase3`.
+        integer :: number_of_constituents = 0
+
+        ! Initialized by `dyn_mpas_define_scalar`.
+        character(strkind), allocatable :: constituent_name(:)
+        integer, allocatable :: index_constituent_to_mpas_scalar(:)
+        integer, allocatable :: index_mpas_scalar_to_constituent(:)
+        logical, allocatable :: is_water_species(:)
+
+        ! Initialized by `dyn_mpas_init_phase4`.
+        integer :: coupling_time_interval
     contains
         private
 
@@ -79,12 +105,24 @@ end subroutine model_error_if
         procedure, pass, public :: read_namelist => dyn_mpas_read_namelist
         procedure, pass, public :: init_phase2 => dyn_mpas_init_phase2
         procedure, pass, public :: init_phase3 => dyn_mpas_init_phase3
+        procedure, pass, public :: define_scalar => dyn_mpas_define_scalar
         procedure, pass, public :: read_write_stream => dyn_mpas_read_write_stream
         procedure, pass :: init_stream_with_pool => dyn_mpas_init_stream_with_pool
+        procedure, pass :: check_variable_status => dyn_mpas_check_variable_status
         procedure, pass, public :: exchange_halo => dyn_mpas_exchange_halo
+        procedure, pass, public :: compute_unit_vector => dyn_mpas_compute_unit_vector
+        procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind
+        procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4
 
         ! Accessor subroutines for users to access 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
+
+        procedure, pass, public :: map_mpas_scalar_index => dyn_mpas_map_mpas_scalar_index
+        procedure, pass, public :: map_constituent_index => dyn_mpas_map_constituent_index
+
+        procedure, pass, public :: get_local_mesh_dimension => dyn_mpas_get_local_mesh_dimension
         procedure, pass, public :: get_global_mesh_dimension => dyn_mpas_get_global_mesh_dimension
 
         procedure, pass :: get_pool_pointer => dyn_mpas_get_pool_pointer
@@ -231,7 +269,6 @@ end subroutine model_error_if
     type(var_info_type), parameter :: input_var_info_list(*) = [ &
         var_info_type('Time'                            , 'real'      , 1), &
         var_info_type('initial_time'                    , 'character' , 0), &
-        var_info_type('relhum'                          , 'real'      , 3), &
         var_info_type('rho'                             , 'real'      , 3), &
         var_info_type('rho_base'                        , 'real'      , 3), &
         var_info_type('scalars'                         , 'real'      , 3), &
@@ -484,7 +521,7 @@ subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_unit, mpas
         ! 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` function pointer.
+        ! 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
@@ -518,12 +555,15 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, &
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_read_namelist'
         character(strkind) :: mpas_calendar
-        character(strkind), pointer :: config_pointer_c => null()
+        character(strkind), pointer :: config_pointer_c
         integer :: ierr
-        logical, pointer :: config_pointer_l => null()
+        logical, pointer :: config_pointer_l
 
         call self % debug_print(subname // ' entered')
 
+        nullify(config_pointer_c)
+        nullify(config_pointer_l)
+
         call self % debug_print('Reading namelist at ', [namelist_path])
 
         ! Override namelist filename so that we can rely on upstream MPAS functionality for reading its own namelist.
@@ -703,37 +743,49 @@ end subroutine dyn_mpas_init_phase2
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_init_phase3
     !
-    !> \brief  Tracks `mpas_init` up to the point of calling `core_init`
+    !> \brief  Tracks `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 point where `core_init` is called. It completes MPAS framework
+    !>  the point where `atm_core_init` is called. It completes MPAS framework
     !>  initialization, including the allocation of all blocks and fields managed
-    !>  by MPAS.
+    !>  by MPAS. However, 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)
     !
     !-------------------------------------------------------------------------------
-    subroutine dyn_mpas_init_phase3(self, cam_pcnst, pio_file)
-        class(mpas_dynamical_core_type), intent(in) :: self
-        integer, intent(in) :: cam_pcnst
+    subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file)
+        class(mpas_dynamical_core_type), intent(inout) :: self
+        integer, intent(in) :: number_of_constituents
         type(file_desc_t), pointer, intent(in) :: pio_file
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase3'
         character(strkind) :: mesh_filename
         integer :: mesh_format
+        integer, pointer :: num_scalars
+        type(mpas_pool_type), pointer :: mpas_pool
 
         call self % debug_print(subname // ' entered')
 
-        call self % debug_print('Number of constituents is ', [cam_pcnst])
+        nullify(mpas_pool)
+        nullify(num_scalars)
+
+        ! In MPAS, there must be at least one constituent, `qv`, which denotes water vapor mixing ratio.
+        ! Because MPAS has some hard-coded array accesses through the `index_qv` index, it will crash
+        ! (i.e., segmentation fault due to invalid memory access) if `qv` is not allocated.
+        self % number_of_constituents = max(1, number_of_constituents)
+
+        call self % debug_print('Number of constituents is ', [self % number_of_constituents])
 
         ! 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`.
-        call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', cam_pcnst)
+        call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', self % number_of_constituents)
 
         ! Not actually used because a PIO file descriptor is directly supplied.
         mesh_filename = 'external mesh'
@@ -759,9 +811,293 @@ subroutine dyn_mpas_init_phase3(self, cam_pcnst, 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.
+        call self % get_pool_pointer(mpas_pool, 'state')
+
+        call mpas_pool_get_dimension(mpas_pool, 'num_scalars', num_scalars)
+
+        if (.not. associated(num_scalars)) then
+            call self % model_error('Failed to find variable "num_scalars"', subname, __LINE__)
+        end if
+
+        ! While we are at it, check if its value is consistent.
+        if (num_scalars /= self % number_of_constituents) then
+            call self % model_error('Failed to allocate constituents', subname, __LINE__)
+        end if
+
+        call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'num_scalars', num_scalars)
+
+        nullify(mpas_pool)
+        nullify(num_scalars)
+
         call self % debug_print(subname // ' completed')
     end subroutine dyn_mpas_init_phase3
 
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_define_scalar
+    !
+    !> \brief  Defines the names of constituents at run-time
+    !> \author Michael Duda
+    !> \date   21 May 2020
+    !> \details
+    !>  Given arrays of constituent names and their corresponding waterness, which
+    !>  must have sizes equal to the number of constituents used to call
+    !>  `dyn_mpas_init_phase3`, this subroutine defines the scalars inside MPAS.
+    !>  Note that MPAS uses the term "scalar", but CAM-SIMA calls it "constituent".
+    !>  Furthermore, because MPAS expects all water scalars to appear in a
+    !>  contiguous index range, this subroutine may reorder the scalars to satisfy
+    !>  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)
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
+        class(mpas_dynamical_core_type), intent(inout) :: self
+        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.
+        character(*), parameter :: mpas_scalar_qv_standard_name(*) = [ character(strkind) :: &
+            'water_vapor_mixing_ratio_wrt_dry_air', &
+            'water_vapor_mixing_ratio_wrt_moist_air', &
+            'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water' &
+        ]
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_define_scalar'
+        integer :: i, j, ierr
+        integer :: index_qv, index_water_start, index_water_end
+        integer :: time_level
+        type(field3dreal), pointer :: field_3d_real
+        type(mpas_pool_type), pointer :: mpas_pool
+
+        call self % debug_print(subname // ' entered')
+
+        nullify(field_3d_real)
+        nullify(mpas_pool)
+
+        if (self % number_of_constituents == 0) then
+            call self % model_error('Constituents must be allocated before being defined', subname, __LINE__)
+        end if
+
+        ! Input sanitization.
+
+        if (size(constituent_name) /= size(is_water_species)) then
+            call self % model_error('Mismatch between numbers of constituent names and their waterness', subname, __LINE__)
+        end if
+
+        if (size(constituent_name) == 0 .and. self % number_of_constituents == 1) then
+            ! If constituent definitions are empty, `qv` is the only constituent per MPAS requirements.
+            ! See `dyn_mpas_init_phase3` for details.
+            allocate(self % constituent_name(1), stat=ierr)
+
+            if (ierr /= 0) then
+                call self % model_error('Failed to allocate constituent_name', subname, __LINE__)
+            end if
+
+            allocate(self % is_water_species(1), stat=ierr)
+
+            if (ierr /= 0) then
+                call self % model_error('Failed to allocate is_water_species', subname, __LINE__)
+            end if
+
+            self % constituent_name(1) = mpas_scalar_qv_standard_name(1)
+            self % is_water_species(1) = .true.
+        else
+            if (size(constituent_name) /= self % number_of_constituents) then
+                call self % model_error('Mismatch between numbers of constituents and their names', subname, __LINE__)
+            end if
+
+            if (any(len_trim(adjustl(constituent_name)) > len(self % constituent_name))) then
+                call self % model_error('Constituent names are too long', subname, __LINE__)
+            end if
+
+            allocate(self % constituent_name(self % number_of_constituents), stat=ierr)
+
+            if (ierr /= 0) then
+                call self % model_error('Failed to allocate constituent_name', subname, __LINE__)
+            end if
+
+            self % constituent_name(:) = adjustl(constituent_name)
+
+            allocate(self % is_water_species(self % number_of_constituents), stat=ierr)
+
+            if (ierr /= 0) then
+                call self % model_error('Failed to allocate is_water_species', subname, __LINE__)
+            end if
+
+            self % is_water_species(:) = is_water_species(:)
+
+            if (size(self % constituent_name) /= size(index_unique(self % constituent_name))) then
+                call self % model_error('Constituent names must be unique', subname, __LINE__)
+            end if
+
+            ! `qv` must be present in constituents per MPAS requirements. It is a water species by definition.
+            ! See `dyn_mpas_init_phase3` for details.
+            index_qv = 0
+
+            ! Lower index in `mpas_scalar_qv_standard_name` has higher precedence, with index 1 being exactly what MPAS wants.
+            set_index_qv: do i = 1, size(mpas_scalar_qv_standard_name)
+                do j = 1, self % number_of_constituents
+                    if (self % constituent_name(j) == mpas_scalar_qv_standard_name(i) .and. self % is_water_species(j)) then
+                        index_qv = j
+
+                        ! The best candidate of `qv` has been found. Exit prematurely.
+                        exit set_index_qv
+                    end if
+                end do
+            end do set_index_qv
+
+            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__)
+            end if
+        end if
+
+        ! Create index mapping between MPAS scalars and constituent names. For example,
+        ! MPAS scalar index `i` corresponds to constituent index `index_mpas_scalar_to_constituent(i)`.
+
+        allocate(self % index_mpas_scalar_to_constituent(self % number_of_constituents), stat=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to allocate index_mpas_scalar_to_constituent', subname, __LINE__)
+        end if
+
+        self % index_mpas_scalar_to_constituent(:) = 0
+        j = 1
+
+        ! Place water species first per MPAS requirements.
+        do i = 1, self % number_of_constituents
+            if (self % is_water_species(i)) then
+                self % index_mpas_scalar_to_constituent(j) = i
+                j = j + 1
+            end if
+        end do
+
+        index_water_start = 1
+        index_water_end = count(self % is_water_species)
+
+        ! Place non-water species second per MPAS requirements.
+        do i = 1, self % number_of_constituents
+            if (.not. self % is_water_species(i)) then
+                self % index_mpas_scalar_to_constituent(j) = i
+                j = j + 1
+            end if
+        end do
+
+        ! Create inverse index mapping between MPAS scalars and constituent names. For example,
+        ! Constituent index `i` corresponds to MPAS scalar index `index_constituent_to_mpas_scalar(i)`.
+
+        allocate(self % index_constituent_to_mpas_scalar(self % number_of_constituents), stat=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to allocate index_constituent_to_mpas_scalar', subname, __LINE__)
+        end if
+
+        self % index_constituent_to_mpas_scalar(:) = 0
+
+        do i = 1, self % number_of_constituents
+            self % index_constituent_to_mpas_scalar(self % index_mpas_scalar_to_constituent(i)) = i
+        end do
+
+        ! Set the index of `qv` in terms of MPAS scalars.
+        index_qv = self % index_constituent_to_mpas_scalar(index_qv)
+
+        ! Print information about constituents.
+        do i = 1, self % number_of_constituents
+            call self % debug_print('Constituent index ' // stringify([i]))
+            call self % debug_print('    Constituent name: ' // &
+                stringify([self % constituent_name(i)]))
+            call self % debug_print('    Is water species: ' // &
+                stringify([self % is_water_species(i)]))
+            call self % debug_print('    Index mapping from constituent to MPAS scalar: ' // &
+                stringify([i]) // ' -> ' // stringify([self % index_constituent_to_mpas_scalar(i)]))
+        end do
+
+        ! Define "scalars" for MPAS.
+
+        call self % get_pool_pointer(mpas_pool, 'state')
+
+        call mpas_pool_add_dimension(mpas_pool, 'index_qv', index_qv)
+        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.
+        time_level = 2
+
+        do i = 1, time_level
+            call mpas_pool_get_field(mpas_pool, 'scalars', field_3d_real, timelevel=i)
+
+            if (.not. associated(field_3d_real)) then
+                call self % model_error('Failed to find variable "scalars"', subname, __LINE__)
+            end if
+
+            do j = 1, self % number_of_constituents
+                field_3d_real % constituentnames(j) = &
+                    trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j))))
+
+                ! Print information about MPAS scalars. Only do it once.
+                if (i == 1) then
+                    call self % debug_print('MPAS scalar index ' // stringify([j]))
+                    call self % debug_print('    MPAS scalar name: ' // &
+                        stringify([field_3d_real % constituentnames(j)]))
+                    call self % debug_print('    Is water species: ' // &
+                        stringify([self % is_water_species(self % index_mpas_scalar_to_constituent(j))]))
+                    call self % debug_print('    Index mapping from MPAS scalar to constituent: ' // &
+                        stringify([j]) // ' -> ' // stringify([self % index_mpas_scalar_to_constituent(j)]))
+                end if
+            end do
+
+            nullify(field_3d_real)
+        end do
+
+        nullify(mpas_pool)
+
+        ! Define "scalars_tend" for MPAS.
+
+        call self % get_pool_pointer(mpas_pool, 'tend')
+
+        call mpas_pool_add_dimension(mpas_pool, 'index_qv', index_qv)
+        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.
+        time_level = 1
+
+        do i = 1, time_level
+            call mpas_pool_get_field(mpas_pool, 'scalars_tend', field_3d_real, timelevel=i)
+
+            if (.not. associated(field_3d_real)) then
+                call self % model_error('Failed to find variable "scalars_tend"', subname, __LINE__)
+            end if
+
+            do j = 1, self % number_of_constituents
+                field_3d_real % constituentnames(j) = &
+                    'tendency_of_' // trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j))))
+            end do
+
+            nullify(field_3d_real)
+        end do
+
+        nullify(mpas_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)
+        call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_end', index_water_end)
+
+        call self % debug_print('index_qv = ' // stringify([index_qv]))
+        call self % debug_print('moist_start = ' // stringify([index_water_start]))
+        call self % debug_print('moist_end = ' // stringify([index_water_end]))
+
+        call self % debug_print(subname // ' completed')
+    end subroutine dyn_mpas_define_scalar
+
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_read_write_stream
     !
@@ -784,12 +1120,15 @@ subroutine dyn_mpas_read_write_stream(self, pio_file, stream_mode, stream_name)
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_read_write_stream'
         integer :: i, ierr
-        type(mpas_pool_type), pointer :: mpas_pool => null()
-        type(mpas_stream_type), pointer :: mpas_stream => null()
+        type(mpas_pool_type), pointer :: mpas_pool
+        type(mpas_stream_type), pointer :: mpas_stream
         type(var_info_type), allocatable :: var_info_list(:)
 
         call self % debug_print(subname // ' entered')
 
+        nullify(mpas_pool)
+        nullify(mpas_stream)
+
         call self % debug_print('Initializing stream "' // trim(adjustl(stream_name)) // '"')
 
         call self % init_stream_with_pool(mpas_pool, mpas_stream, pio_file, stream_mode, stream_name)
@@ -894,23 +1233,40 @@ 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, ndims, stream_format, varid
-        type(field0dchar), pointer :: field_0d_char => null()
-        type(field1dchar), pointer :: field_1d_char => null()
-        type(field0dinteger), pointer :: field_0d_integer => null()
-        type(field1dinteger), pointer :: field_1d_integer => null()
-        type(field2dinteger), pointer :: field_2d_integer => null()
-        type(field3dinteger), pointer :: field_3d_integer => null()
-        type(field0dreal), pointer :: field_0d_real => null()
-        type(field1dreal), pointer :: field_1d_real => null()
-        type(field2dreal), pointer :: field_2d_real => null()
-        type(field3dreal), pointer :: field_3d_real => null()
-        type(field4dreal), pointer :: field_4d_real => null()
-        type(field5dreal), pointer :: field_5d_real => null()
+        integer :: i, ierr, stream_format
+        ! 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`).
+        logical, allocatable :: var_is_tkr_compatible(:)
+        type(field0dchar), pointer :: field_0d_char
+        type(field1dchar), pointer :: field_1d_char
+        type(field0dinteger), pointer :: field_0d_integer
+        type(field1dinteger), pointer :: field_1d_integer
+        type(field2dinteger), pointer :: field_2d_integer
+        type(field3dinteger), pointer :: field_3d_integer
+        type(field0dreal), pointer :: field_0d_real
+        type(field1dreal), pointer :: field_1d_real
+        type(field2dreal), pointer :: field_2d_real
+        type(field3dreal), pointer :: field_3d_real
+        type(field4dreal), pointer :: field_4d_real
+        type(field5dreal), pointer :: field_5d_real
         type(var_info_type), allocatable :: var_info_list(:)
 
         call self % debug_print(subname // ' entered')
 
+        nullify(field_0d_char)
+        nullify(field_1d_char)
+        nullify(field_0d_integer)
+        nullify(field_1d_integer)
+        nullify(field_2d_integer)
+        nullify(field_3d_integer)
+        nullify(field_0d_real)
+        nullify(field_1d_real)
+        nullify(field_2d_real)
+        nullify(field_3d_real)
+        nullify(field_4d_real)
+        nullify(field_5d_real)
+
         call mpas_pool_create_pool(mpas_pool)
 
         allocate(mpas_stream, stat=ierr)
@@ -970,28 +1326,23 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
                 stringify([var_info_list(i) % rank]))
 
             if (trim(adjustl(stream_mode)) == 'r' .or. trim(adjustl(stream_mode)) == 'read') then
-                ! Check if "<variable name>" is present.
-                ierr = pio_inq_varid(pio_file, trim(adjustl(var_info_list(i) % name)), varid)
+                call self % check_variable_status(var_is_present, var_is_tkr_compatible, pio_file, var_info_list(i))
 
                 ! Do not hard crash the model if a variable is missing and cannot be read.
                 ! This can happen if users attempt to initialize/restart the model with data generated by
                 ! older versions of MPAS. Print a debug message to let users decide if this is acceptable.
-                if (ierr /= pio_noerr) then
-                    call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // '"')
+                if (.not. any(var_is_present)) then
+                    call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // &
+                        '" due to not present')
 
                     cycle
                 end if
 
-                ierr = pio_inq_varndims(pio_file, varid, ndims)
-
-                if (ierr /= pio_noerr) then
-                    call self % model_error('Failed to inquire variable rank for "' // trim(adjustl(var_info_list(i) % name)) // &
-                        '"', subname, __LINE__)
-                end if
+                if (any(var_is_present .and. .not. var_is_tkr_compatible)) then
+                    call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // &
+                        '" due to not TKR compatible')
 
-                if (ndims /= var_info_list(i) % rank) then
-                    call self % model_error('Variable rank mismatch for "' // trim(adjustl(var_info_list(i) % name)) // &
-                        '"', subname, __LINE__)
+                    cycle
                 end if
             end if
 
@@ -1031,8 +1382,8 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
 
                             nullify(field_1d_char)
                         case default
-                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), &
-                                subname, __LINE__)
+                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // &
+                                ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__)
                     end select
                 case ('integer')
                     select case (var_info_list(i) % rank)
@@ -1085,8 +1436,8 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
 
                             nullify(field_3d_integer)
                         case default
-                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), &
-                                subname, __LINE__)
+                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // &
+                                ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__)
                     end select
                 case ('real')
                     select case (var_info_list(i) % rank)
@@ -1163,12 +1514,12 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
 
                             nullify(field_5d_real)
                         case default
-                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), &
-                                subname, __LINE__)
+                            call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // &
+                                ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__)
                     end select
                 case default
                     call self % model_error('Unsupported variable type "' // trim(adjustl(var_info_list(i) % type)) // &
-                        '"', subname, __LINE__)
+                        '" for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__)
             end select
 
             if (ierr /= mpas_stream_noerr) then
@@ -1177,8 +1528,6 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
             end if
         end do
 
-        deallocate(var_info_list)
-
         if (trim(adjustl(stream_mode)) == 'w' .or. trim(adjustl(stream_mode)) == 'write') then
             ! Add MPAS-specific attributes to stream.
             call self % debug_print('Adding attributes to stream')
@@ -1194,7 +1543,7 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file
             call add_stream_attribute('is_periodic', self % domain_ptr % is_periodic)
             call add_stream_attribute('mesh_spec', self % domain_ptr % mesh_spec)
             call add_stream_attribute('on_a_sphere', self % domain_ptr % on_a_sphere)
-            call add_stream_attribute('parent_id', self % domain_ptr %  parent_id)
+            call add_stream_attribute('parent_id', self % domain_ptr % parent_id)
             call add_stream_attribute('sphere_radius', self % domain_ptr % sphere_radius)
             call add_stream_attribute('x_period', self % domain_ptr % x_period)
             call add_stream_attribute('y_period', self % domain_ptr % y_period)
@@ -1264,70 +1613,141 @@ subroutine add_stream_attribute_1d(attribute_name, attribute_value)
         end subroutine add_stream_attribute_1d
     end subroutine dyn_mpas_init_stream_with_pool
 
-    !> Parse one or more stream names and return variable information contained in those streams as a list of `var_info_type`.
-    !> Multiple stream names should be separated by `+` (i.e., a plus).
+    !> 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.
-    !> (KCW, 2024-03-15)
-    pure recursive function parse_stream_name(stream_name) result(var_info_list)
+    !> (KCW, 2024-06-01)
+    pure function parse_stream_name(stream_name) result(var_info_list)
         character(*), intent(in) :: stream_name
         type(var_info_type), allocatable :: var_info_list(:)
 
-        character(64), allocatable :: var_name_list(:)
-        integer :: i, n, offset
-        type(var_info_type), allocatable :: var_info_list_append(:)
+        character(*), parameter :: supported_stream_name_operator = '+-'
+        character(1) :: stream_name_operator
+        character(:), allocatable :: stream_name_fragment
+        character(len(invariant_var_info_list % name)), allocatable :: var_name_list(:)
+        integer :: i, j, n, offset
+        type(var_info_type), allocatable :: var_info_list_buffer(:)
 
-        allocate(var_info_list(0))
+        n = len_trim(stream_name)
 
-        n = len(stream_name)
-        offset = 0
+        if (n == 0) then
+            ! Empty character string means empty list.
+            var_info_list = parse_stream_name_fragment('')
+
+            return
+        end if
+
+        i = scan(stream_name, supported_stream_name_operator)
+
+        if (i == 0) then
+            ! No operators are present in the stream name. It is just a single stream name fragment.
+            stream_name_fragment = stream_name
+            var_info_list = parse_stream_name_fragment(stream_name_fragment)
 
-        if (offset + 1 > n) then
             return
         end if
 
-        i = index(stream_name(offset + 1:), '+')
+        offset = 0
+        var_info_list = parse_stream_name_fragment('')
+
+        do while (.true.)
+            ! Extract operator from the stream name.
+            if (offset > 0) then
+                stream_name_operator = stream_name(offset:offset)
+            else
+                stream_name_operator = '+'
+            end if
 
-        do while (i > 0)
+            ! Extract stream name fragment from the stream name.
             if (i > 1) then
-                var_info_list_append = parse_stream_name(stream_name(offset + 1:offset + i - 1))
-                var_info_list = [var_info_list, var_info_list_append]
+                stream_name_fragment = stream_name(offset + 1:offset + i - 1)
+            else
+                stream_name_fragment = ''
+            end if
 
-                deallocate(var_info_list_append)
+            ! Process the stream name fragment according to the operator.
+            if (len_trim(stream_name_fragment) > 0) then
+                var_info_list_buffer = parse_stream_name_fragment(stream_name_fragment)
+
+                select case (stream_name_operator)
+                    case ('+')
+                        var_info_list = [var_info_list, var_info_list_buffer]
+                    case ('-')
+                        do j = 1, size(var_info_list_buffer)
+                            var_name_list = var_info_list % name
+                            var_info_list = pack(var_info_list, var_name_list /= var_info_list_buffer(j) % name)
+                        end do
+                    case default
+                        ! Do nothing for unknown operators. Should not happen at all.
+                end select
             end if
 
             offset = offset + i
 
+            ! Terminate loop when everything in the stream name has been processed.
             if (offset + 1 > n) then
                 exit
             end if
 
-            i = index(stream_name(offset + 1:), '+')
+            i = scan(stream_name(offset + 1:), supported_stream_name_operator)
+
+            ! Run the loop one last time for the remaining stream name fragment.
+            if (i == 0) then
+                i = n - offset + 1
+            end if
         end do
 
-        if (offset + 1 > n) then
-            return
-        end if
+        ! Discard duplicate variable information by names.
+        var_name_list = var_info_list % name
+        var_info_list = var_info_list(index_unique(var_name_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.
+    !> (KCW, 2024-06-01)
+    pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_list)
+        character(*), intent(in) :: stream_name_fragment
+        type(var_info_type), allocatable :: var_info_list(:)
 
-        select case (trim(adjustl(stream_name(offset + 1:))))
+        character(len(invariant_var_info_list % name)), allocatable :: var_name_list(:)
+        type(var_info_type), allocatable :: var_info_list_buffer(:)
+
+        select case (trim(adjustl(stream_name_fragment)))
+            case ('')
+                allocate(var_info_list(0))
             case ('invariant')
-                allocate(var_info_list_append, source=invariant_var_info_list)
+                allocate(var_info_list, source=invariant_var_info_list)
             case ('input')
-                allocate(var_info_list_append, source=input_var_info_list)
+                allocate(var_info_list, source=input_var_info_list)
             case ('restart')
-                allocate(var_info_list_append, source=restart_var_info_list)
+                allocate(var_info_list, source=restart_var_info_list)
             case default
-                allocate(var_info_list_append(0))
-        end select
+                allocate(var_info_list(0))
 
-        var_info_list = [var_info_list, var_info_list_append]
+                var_name_list = invariant_var_info_list % name
 
-        ! Discard duplicate variable information by names.
-        var_name_list = var_info_list(:) % name
-        var_info_list = var_info_list(index_unique(var_name_list))
+                if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then
+                    var_info_list_buffer = pack(invariant_var_info_list, var_name_list == trim(adjustl(stream_name_fragment)))
+                    var_info_list = [var_info_list, var_info_list_buffer]
+                end if
 
-        deallocate(var_info_list_append)
-        deallocate(var_name_list)
-    end function parse_stream_name
+                var_name_list = input_var_info_list % name
+
+                if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then
+                    var_info_list_buffer = pack(input_var_info_list, var_name_list == trim(adjustl(stream_name_fragment)))
+                    var_info_list = [var_info_list, var_info_list_buffer]
+                end if
+
+                var_name_list = restart_var_info_list % name
+
+                if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then
+                    var_info_list_buffer = pack(restart_var_info_list, var_name_list == trim(adjustl(stream_name_fragment)))
+                    var_info_list = [var_info_list, var_info_list_buffer]
+                end if
+        end select
+    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.
@@ -1408,170 +1828,1220 @@ pure function index_unique(array)
     end function index_unique
 
     !-------------------------------------------------------------------------------
-    ! subroutine dyn_mpas_exchange_halo
+    ! subroutine dyn_mpas_check_variable_status
     !
-    !> \brief  Updates the halo layers of the named field
-    !> \author Michael Duda
-    !> \date   16 January 2020
+    !> \brief  Check and return variable status on the given file
+    !> \author Kuan-Chih Wang
+    !> \date   2024-06-04
     !> \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)
+    !>  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.
+    !>  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
+    !>  itself. Otherwise, for a variable array, the checks are performed on its
+    !>  constituent parts instead.
     !
     !-------------------------------------------------------------------------------
-    subroutine dyn_mpas_exchange_halo(self, field_name)
+    subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compatible, pio_file, var_info)
         class(mpas_dynamical_core_type), intent(in) :: self
-        character(*), intent(in) :: field_name
-
-        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_exchange_halo'
-        type(field1dinteger), pointer :: field_1d_integer => null()
-        type(field2dinteger), pointer :: field_2d_integer => null()
-        type(field3dinteger), pointer :: field_3d_integer => null()
-        type(field1dreal), pointer :: field_1d_real => null()
-        type(field2dreal), pointer :: field_2d_real => null()
-        type(field3dreal), pointer :: field_3d_real => null()
-        type(field4dreal), pointer :: field_4d_real => null()
-        type(field5dreal), pointer :: field_5d_real => null()
-        type(mpas_pool_field_info_type) :: mpas_pool_field_info
+        logical, allocatable, intent(out) :: var_is_present(:)
+        logical, allocatable, intent(out) :: var_is_tkr_compatible(:)
+        type(file_desc_t), pointer, intent(in) :: pio_file
+        type(var_info_type), intent(in) :: var_info
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_check_variable_status'
+        character(strkind), allocatable :: var_name_list(:)
+        integer :: i, ierr, varid, varndims, vartype
+        type(field0dchar), pointer :: field_0d_char
+        type(field1dchar), pointer :: field_1d_char
+        type(field0dinteger), pointer :: field_0d_integer
+        type(field1dinteger), pointer :: field_1d_integer
+        type(field2dinteger), pointer :: field_2d_integer
+        type(field3dinteger), pointer :: field_3d_integer
+        type(field0dreal), pointer :: field_0d_real
+        type(field1dreal), pointer :: field_1d_real
+        type(field2dreal), pointer :: field_2d_real
+        type(field3dreal), pointer :: field_3d_real
+        type(field4dreal), pointer :: field_4d_real
+        type(field5dreal), pointer :: field_5d_real
 
         call self % debug_print(subname // ' entered')
 
-        call self % debug_print('Inquiring field information for "' // trim(adjustl(field_name)) // '"')
+        nullify(field_0d_char)
+        nullify(field_1d_char)
+        nullify(field_0d_integer)
+        nullify(field_1d_integer)
+        nullify(field_2d_integer)
+        nullify(field_3d_integer)
+        nullify(field_0d_real)
+        nullify(field_1d_real)
+        nullify(field_2d_real)
+        nullify(field_3d_real)
+        nullify(field_4d_real)
+        nullify(field_5d_real)
+
+        ! Extract a list of variable names to check on the file.
+        ! For an ordinary variable, this list just contains its name.
+        ! For a variable array, this list contains the names of its constituent parts.
+        select case (trim(adjustl(var_info % type)))
+            case ('character')
+                select case (var_info % rank)
+                    case (0)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_0d_char, timelevel=1)
 
-        call mpas_pool_get_field_info(self % domain_ptr % blocklist % allfields, &
-            trim(adjustl(field_name)), mpas_pool_field_info)
+                        if (.not. associated(field_0d_char)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
 
-        if (mpas_pool_field_info % fieldtype == -1 .or. &
-            mpas_pool_field_info % ndims == -1 .or. &
-            mpas_pool_field_info % nhalolayers == -1) then
-            call self % model_error('Invalid field information for "' // trim(adjustl(field_name)) // '"', subname, __LINE__)
-        end if
+                        if (field_0d_char % isvararray .and. associated(field_0d_char % constituentnames)) then
+                            allocate(var_name_list(size(field_0d_char % constituentnames)), stat=ierr)
 
-        ! No halo layers to exchange. This field is not decomposed.
-        if (mpas_pool_field_info % nhalolayers == 0) then
-            call self % debug_print('Skipping field "' // trim(adjustl(field_name)) // '"')
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
 
-            return
-        end if
+                            var_name_list(:) = field_0d_char % constituentnames(:)
+                        end if
 
-        call self % debug_print('Exchanging halo layers for "' // trim(adjustl(field_name)) // '"')
+                        nullify(field_0d_char)
+                    case (1)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_1d_char, timelevel=1)
 
-        select case (mpas_pool_field_info % fieldtype)
-            case (mpas_pool_integer)
-                select case (mpas_pool_field_info % ndims)
+                        if (.not. associated(field_1d_char)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_1d_char % isvararray .and. associated(field_1d_char % constituentnames)) then
+                            allocate(var_name_list(size(field_1d_char % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_1d_char % constituentnames(:)
+                        end if
+
+                        nullify(field_1d_char)
+                    case default
+                        call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // &
+                            ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__)
+                end select
+            case ('integer')
+                select case (var_info % rank)
+                    case (0)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_0d_integer, timelevel=1)
+
+                        if (.not. associated(field_0d_integer)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_0d_integer % isvararray .and. associated(field_0d_integer % constituentnames)) then
+                            allocate(var_name_list(size(field_0d_integer % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_0d_integer % constituentnames(:)
+                        end if
+
+                        nullify(field_0d_integer)
                     case (1)
                         call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_1d_integer, timelevel=1)
+                            trim(adjustl(var_info % name)), field_1d_integer, timelevel=1)
 
                         if (.not. associated(field_1d_integer)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
                                 '"', subname, __LINE__)
                         end if
 
-                        call mpas_dmpar_exch_halo_field(field_1d_integer)
+                        if (field_1d_integer % isvararray .and. associated(field_1d_integer % constituentnames)) then
+                            allocate(var_name_list(size(field_1d_integer % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_1d_integer % constituentnames(:)
+                        end if
 
                         nullify(field_1d_integer)
                     case (2)
                         call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_2d_integer, timelevel=1)
+                            trim(adjustl(var_info % name)), field_2d_integer, timelevel=1)
 
                         if (.not. associated(field_2d_integer)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
                                 '"', subname, __LINE__)
                         end if
 
-                        call mpas_dmpar_exch_halo_field(field_2d_integer)
+                        if (field_2d_integer % isvararray .and. associated(field_2d_integer % constituentnames)) then
+                            allocate(var_name_list(size(field_2d_integer % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_2d_integer % constituentnames(:)
+                        end if
 
                         nullify(field_2d_integer)
                     case (3)
                         call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_3d_integer, timelevel=1)
+                            trim(adjustl(var_info % name)), field_3d_integer, timelevel=1)
 
                         if (.not. associated(field_3d_integer)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
                                 '"', subname, __LINE__)
                         end if
 
-                        call mpas_dmpar_exch_halo_field(field_3d_integer)
+                        if (field_3d_integer % isvararray .and. associated(field_3d_integer % constituentnames)) then
+                            allocate(var_name_list(size(field_3d_integer % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_3d_integer % constituentnames(:)
+                        end if
 
                         nullify(field_3d_integer)
                     case default
-                        call self % model_error('Unsupported field rank ' // stringify([mpas_pool_field_info % ndims]), &
-                            subname, __LINE__)
+                        call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // &
+                            ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__)
                 end select
-            case (mpas_pool_real)
-                select case (mpas_pool_field_info % ndims)
-                    case (1)
+            case ('real')
+                select case (var_info % rank)
+                    case (0)
                         call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_1d_real, timelevel=1)
+                            trim(adjustl(var_info % name)), field_0d_real, timelevel=1)
 
-                        if (.not. associated(field_1d_real)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                        if (.not. associated(field_0d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
                                 '"', subname, __LINE__)
                         end if
 
-                        call mpas_dmpar_exch_halo_field(field_1d_real)
+                        if (field_0d_real % isvararray .and. associated(field_0d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_0d_real % constituentnames)), stat=ierr)
 
-                        nullify(field_1d_real)
-                    case (2)
-                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_0d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_0d_real)
+                    case (1)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_1d_real, timelevel=1)
+
+                        if (.not. associated(field_1d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_1d_real % isvararray .and. associated(field_1d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_1d_real % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_1d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_1d_real)
+                    case (2)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_2d_real, timelevel=1)
+
+                        if (.not. associated(field_2d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_2d_real % isvararray .and. associated(field_2d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_2d_real % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_2d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_2d_real)
+                    case (3)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_3d_real, timelevel=1)
+
+                        if (.not. associated(field_3d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_3d_real % isvararray .and. associated(field_3d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_3d_real % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_3d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_3d_real)
+                    case (4)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_4d_real, timelevel=1)
+
+                        if (.not. associated(field_4d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_4d_real % isvararray .and. associated(field_4d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_4d_real % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_4d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_4d_real)
+                    case (5)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(var_info % name)), field_5d_real, timelevel=1)
+
+                        if (.not. associated(field_5d_real)) then
+                            call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        if (field_5d_real % isvararray .and. associated(field_5d_real % constituentnames)) then
+                            allocate(var_name_list(size(field_5d_real % constituentnames)), stat=ierr)
+
+                            if (ierr /= 0) then
+                                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+                            end if
+
+                            var_name_list(:) = field_5d_real % constituentnames(:)
+                        end if
+
+                        nullify(field_5d_real)
+                    case default
+                        call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // &
+                            ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__)
+                end select
+            case default
+                call self % model_error('Unsupported variable type "' // trim(adjustl(var_info % type)) // &
+                    '" for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__)
+        end select
+
+        if (.not. allocated(var_name_list)) then
+            allocate(var_name_list(1), stat=ierr)
+
+            if (ierr /= 0) then
+                call self % model_error('Failed to allocate var_name_list', subname, __LINE__)
+            end if
+
+            var_name_list(1) = var_info % name
+        end if
+
+        allocate(var_is_present(size(var_name_list)), stat=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to allocate var_is_present', subname, __LINE__)
+        end if
+
+        var_is_present(:) = .false.
+
+        allocate(var_is_tkr_compatible(size(var_name_list)), stat=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to allocate var_is_tkr_compatible', subname, __LINE__)
+        end if
+
+        var_is_tkr_compatible(:) = .false.
+
+        if (.not. associated(pio_file)) then
+            return
+        end if
+
+        if (.not. pio_file_is_open(pio_file)) then
+            return
+        end if
+
+        do i = 1, size(var_name_list)
+            ! Check if the variable is present on the file.
+            ierr = pio_inq_varid(pio_file, trim(adjustl(var_name_list(i))), varid)
+
+            if (ierr /= pio_noerr) then
+                cycle
+            end if
+
+            var_is_present(i) = .true.
+
+            ! Check if the variable is "TK"R compatible between MPAS and the file.
+            ierr = pio_inq_vartype(pio_file, varid, vartype)
+
+            if (ierr /= pio_noerr) then
+                cycle
+            end if
+
+            select case (trim(adjustl(var_info % type)))
+                case ('character')
+                    if (vartype /= pio_char) then
+                        cycle
+                    end if
+                case ('integer')
+                    if (vartype /= pio_int) then
+                        cycle
+                    end if
+                case ('real')
+                    if (rkind == r4kind .and. vartype /= pio_real) then
+                        cycle
+                    end if
+
+                    if (rkind == r8kind .and. vartype /= pio_double) then
+                        cycle
+                    end if
+                case default
+                    cycle
+            end select
+
+            ! Check if the variable is TK"R" compatible between MPAS and the file.
+            ierr = pio_inq_varndims(pio_file, varid, varndims)
+
+            if (ierr /= pio_noerr) then
+                cycle
+            end if
+
+            if (varndims /= var_info % rank) then
+                cycle
+            end if
+
+            var_is_tkr_compatible(i) = .true.
+        end do
+
+        call self % debug_print('var_name_list = ' // stringify(var_name_list))
+        call self % debug_print('var_is_present = ' // stringify(var_is_present))
+        call self % debug_print('var_is_tkr_compatible = ' // stringify(var_is_tkr_compatible))
+
+        call self % debug_print(subname // ' completed')
+    end subroutine dyn_mpas_check_variable_status
+
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_exchange_halo
+    !
+    !> \brief  Updates 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)
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_exchange_halo(self, field_name)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        character(*), intent(in) :: field_name
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_exchange_halo'
+        type(field1dinteger), pointer :: field_1d_integer
+        type(field2dinteger), pointer :: field_2d_integer
+        type(field3dinteger), pointer :: field_3d_integer
+        type(field1dreal), pointer :: field_1d_real
+        type(field2dreal), pointer :: field_2d_real
+        type(field3dreal), pointer :: field_3d_real
+        type(field4dreal), pointer :: field_4d_real
+        type(field5dreal), pointer :: field_5d_real
+        type(mpas_pool_field_info_type) :: mpas_pool_field_info
+
+        call self % debug_print(subname // ' entered')
+
+        nullify(field_1d_integer)
+        nullify(field_2d_integer)
+        nullify(field_3d_integer)
+        nullify(field_1d_real)
+        nullify(field_2d_real)
+        nullify(field_3d_real)
+        nullify(field_4d_real)
+        nullify(field_5d_real)
+
+        call self % debug_print('Inquiring field information for "' // trim(adjustl(field_name)) // '"')
+
+        call mpas_pool_get_field_info(self % domain_ptr % blocklist % allfields, &
+            trim(adjustl(field_name)), mpas_pool_field_info)
+
+        if (mpas_pool_field_info % fieldtype == -1 .or. &
+            mpas_pool_field_info % ndims == -1 .or. &
+            mpas_pool_field_info % nhalolayers == -1) then
+            call self % model_error('Invalid field information for "' // trim(adjustl(field_name)) // '"', subname, __LINE__)
+        end if
+
+        ! No halo layers to exchange. This field is not decomposed.
+        if (mpas_pool_field_info % nhalolayers == 0) then
+            call self % debug_print('Skipping field "' // trim(adjustl(field_name)) // '"')
+
+            return
+        end if
+
+        call self % debug_print('Exchanging halo layers for "' // trim(adjustl(field_name)) // '"')
+
+        select case (mpas_pool_field_info % fieldtype)
+            case (mpas_pool_integer)
+                select case (mpas_pool_field_info % ndims)
+                    case (1)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_1d_integer, timelevel=1)
+
+                        if (.not. associated(field_1d_integer)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_1d_integer)
+
+                        nullify(field_1d_integer)
+                    case (2)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_2d_integer, timelevel=1)
+
+                        if (.not. associated(field_2d_integer)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_2d_integer)
+
+                        nullify(field_2d_integer)
+                    case (3)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_3d_integer, timelevel=1)
+
+                        if (.not. associated(field_3d_integer)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_3d_integer)
+
+                        nullify(field_3d_integer)
+                    case default
+                        call self % model_error('Unsupported field rank ' // stringify([mpas_pool_field_info % ndims]), &
+                            subname, __LINE__)
+                end select
+            case (mpas_pool_real)
+                select case (mpas_pool_field_info % ndims)
+                    case (1)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_1d_real, timelevel=1)
+
+                        if (.not. associated(field_1d_real)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_1d_real)
+
+                        nullify(field_1d_real)
+                    case (2)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
                             trim(adjustl(field_name)), field_2d_real, timelevel=1)
 
-                        if (.not. associated(field_2d_real)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
-                                '"', subname, __LINE__)
-                        end if
+                        if (.not. associated(field_2d_real)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_2d_real)
+
+                        nullify(field_2d_real)
+                    case (3)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_3d_real, timelevel=1)
+
+                        if (.not. associated(field_3d_real)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_3d_real)
+
+                        nullify(field_3d_real)
+                    case (4)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_4d_real, timelevel=1)
+
+                        if (.not. associated(field_4d_real)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_4d_real)
+
+                        nullify(field_4d_real)
+                    case (5)
+                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
+                            trim(adjustl(field_name)), field_5d_real, timelevel=1)
+
+                        if (.not. associated(field_5d_real)) then
+                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
+                                '"', subname, __LINE__)
+                        end if
+
+                        call mpas_dmpar_exch_halo_field(field_5d_real)
+
+                        nullify(field_5d_real)
+                    case default
+                        call self % model_error('Unsupported field rank ' // stringify([mpas_pool_field_info % ndims]), &
+                            subname, __LINE__)
+                end select
+            case default
+                call self % model_error('Unsupported field type (Must be one of: integer, real)', subname, __LINE__)
+        end select
+
+        call self % debug_print(subname // ' completed')
+    end subroutine dyn_mpas_exchange_halo
+
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_compute_unit_vector
+    !
+    !> \brief  Computes 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
+    !>  `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
+    !>  Ported and refactored for CAM-SIMA. (KCW, 2024-04-23)
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_compute_unit_vector(self)
+        class(mpas_dynamical_core_type), intent(in) :: self
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector'
+        integer :: i
+        integer, pointer :: ncells
+        real(rkind), pointer :: latcell(:), loncell(:)
+        real(rkind), pointer :: east(:, :), north(:, :)
+        type(mpas_pool_type), pointer :: mpas_pool
+
+        call self % debug_print(subname // ' entered')
+
+        nullify(ncells)
+        nullify(latcell, loncell)
+
+        nullify(east, north)
+
+        nullify(mpas_pool)
+
+        ! Input.
+        call self % get_variable_pointer(ncells, 'dim', 'nCells')
+        call self % get_variable_pointer(latcell, 'mesh', 'latCell')
+        call self % get_variable_pointer(loncell, 'mesh', 'lonCell')
+
+        ! Output.
+        call self % get_variable_pointer(east, 'mesh', 'east')
+        call self % get_variable_pointer(north, 'mesh', 'north')
+
+        do i = 1, ncells
+            east(1, i) = -sin(loncell(i))
+            east(2, i) =  cos(loncell(i))
+            east(3, i) =  0.0_rkind
+            ! `r3_normalize` has been inlined below.
+            east(1:3, i) = east(1:3, i) / sqrt(sum(east(1:3, i) * east(1:3, i)))
+
+            north(1, i) = -cos(loncell(i)) * sin(latcell(i))
+            north(2, i) = -sin(loncell(i)) * sin(latcell(i))
+            north(3, i) =  cos(latcell(i))
+            ! `r3_normalize` has been inlined below.
+            north(1:3, i) = north(1:3, i) / sqrt(sum(north(1:3, i) * north(1:3, i)))
+        end do
+
+        nullify(ncells)
+        nullify(latcell, loncell)
+
+        nullify(east, north)
+
+        call self % get_pool_pointer(mpas_pool, 'mesh')
+        call mpas_initialize_vectors(mpas_pool)
+
+        nullify(mpas_pool)
+
+        call self % debug_print(subname // ' completed')
+    end subroutine dyn_mpas_compute_unit_vector
+
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_compute_edge_wind
+    !
+    !> \brief  Computes the edge-normal wind 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 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.
+    !> \addenda
+    !>  Ported and refactored for CAM-SIMA. (KCW, 2024-05-08)
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_compute_edge_wind(self)
+        class(mpas_dynamical_core_type), intent(in) :: self
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_edge_wind'
+        integer :: cell1, cell2, i
+        integer, pointer :: cellsonedge(:, :)
+        integer, pointer :: nedges
+        real(rkind), pointer :: east(:, :), north(:, :)
+        real(rkind), pointer :: edgenormalvectors(:, :)
+        real(rkind), pointer :: ucellzonal(:, :), ucellmeridional(:, :)
+        real(rkind), pointer :: uedge(:, :)
+
+        call self % debug_print(subname // ' entered')
+
+        nullify(nedges)
+
+        nullify(ucellzonal, ucellmeridional)
+
+        nullify(cellsonedge)
+        nullify(east, north)
+        nullify(edgenormalvectors)
+
+        nullify(uedge)
+
+        ! Make sure halo layers are up-to-date before computation.
+        call self % exchange_halo('uReconstructZonal')
+        call self % exchange_halo('uReconstructMeridional')
+
+        ! Input.
+        call self % get_variable_pointer(nedges, 'dim', 'nEdges')
+
+        call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
+        call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
+
+        call self % get_variable_pointer(cellsonedge, 'mesh', 'cellsOnEdge')
+        call self % get_variable_pointer(east, 'mesh', 'east')
+        call self % get_variable_pointer(north, 'mesh', 'north')
+        call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors')
+
+        ! Output.
+        call self % get_variable_pointer(uedge, 'state', 'u', time_level=1)
+
+        do i = 1, nedges
+            cell1 = cellsonedge(1, i)
+            cell2 = cellsonedge(2, i)
+
+            uedge(:, i) = ucellzonal(:, cell1) * 0.5_rkind * (edgenormalvectors(1, i) * east(1, cell1)  + &
+                                                              edgenormalvectors(2, i) * east(2, cell1)  + &
+                                                              edgenormalvectors(3, i) * east(3, cell1)) + &
+                          ucellmeridional(:, cell1) * 0.5_rkind * (edgenormalvectors(1, i) * north(1, cell1)  + &
+                                                                   edgenormalvectors(2, i) * north(2, cell1)  + &
+                                                                   edgenormalvectors(3, i) * north(3, cell1)) + &
+                          ucellzonal(:, cell2) * 0.5_rkind * (edgenormalvectors(1, i) * east(1, cell2)  + &
+                                                              edgenormalvectors(2, i) * east(2, cell2)  + &
+                                                              edgenormalvectors(3, i) * east(3, cell2)) + &
+                          ucellmeridional(:, cell2) * 0.5_rkind * (edgenormalvectors(1, i) * north(1, cell2)  + &
+                                                                   edgenormalvectors(2, i) * north(2, cell2)  + &
+                                                                   edgenormalvectors(3, i) * north(3, cell2))
+        end do
+
+        nullify(nedges)
+
+        nullify(ucellzonal, ucellmeridional)
+
+        nullify(cellsonedge)
+        nullify(east, north)
+        nullify(edgenormalvectors)
+
+        nullify(uedge)
+
+        ! Make sure halo layers are up-to-date after computation.
+        call self % exchange_halo('u')
+
+        call self % debug_print(subname // ' completed')
+    end subroutine dyn_mpas_compute_edge_wind
+
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_init_phase4
+    !
+    !> \brief  Tracks `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)
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
+        class(mpas_dynamical_core_type), intent(inout) :: self
+        integer, intent(in) :: coupling_time_interval ! Set the time interval, in seconds, over which MPAS dynamical core
+                                                      ! should integrate each time it is called to run.
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase4'
+        character(strkind) :: date_time
+        character(strkind), pointer :: initial_time_1, initial_time_2
+        character(strkind), pointer :: xtime
+        integer :: ierr
+        integer, pointer :: nvertlevels, maxedges, maxedges2, num_scalars
+        logical, pointer :: config_do_restart
+        real(rkind), pointer :: config_dt
+        type(field0dreal), pointer :: field_0d_real
+        type(mpas_pool_type), pointer :: mpas_pool
+        type(mpas_time_type) :: mpas_time
+
+        call self % debug_print(subname // ' entered')
+
+        nullify(initial_time_1, initial_time_2)
+        nullify(xtime)
+        nullify(nvertlevels, maxedges, maxedges2, num_scalars)
+        nullify(config_do_restart)
+        nullify(config_dt)
+        nullify(field_0d_real)
+        nullify(mpas_pool)
+
+        if (coupling_time_interval <= 0) then
+            call self % model_error('Invalid coupling time interval ' // stringify([real(coupling_time_interval, rkind)]), &
+                subname, __LINE__)
+        end if
+
+        call self % get_variable_pointer(config_dt, 'cfg', 'config_dt')
 
-                        call mpas_dmpar_exch_halo_field(field_2d_real)
+        if (config_dt <= 0.0_rkind) then
+            call self % model_error('Invalid time step ' // stringify([config_dt]), &
+                subname, __LINE__)
+        end if
 
-                        nullify(field_2d_real)
-                    case (3)
-                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_3d_real, timelevel=1)
+        ! `config_dt` in MPAS is a floating-point number. Testing floating-point numbers for divisibility is not trivial and
+        ! should be done carefully.
+        if (.not. almost_divisible(real(coupling_time_interval, rkind), config_dt)) then
+            call self % model_error('Coupling time interval ' // stringify([real(coupling_time_interval, rkind)]) // &
+                ' must be divisible by time step ' // stringify([config_dt]), subname, __LINE__)
+        end if
 
-                        if (.not. associated(field_3d_real)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
-                                '"', subname, __LINE__)
-                        end if
+        self % coupling_time_interval = coupling_time_interval
 
-                        call mpas_dmpar_exch_halo_field(field_3d_real)
+        call self % debug_print('Coupling time interval is ' // stringify([real(self % coupling_time_interval, rkind)]) // &
+            ' seconds')
+        call self % debug_print('Time step is ' // stringify([config_dt]) // ' seconds')
 
-                        nullify(field_3d_real)
-                    case (4)
-                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_4d_real, timelevel=1)
+        nullify(config_dt)
 
-                        if (.not. associated(field_4d_real)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
-                                '"', subname, __LINE__)
-                        end if
+        ! Compute derived constants.
+        call mpas_constants_compute_derived()
 
-                        call mpas_dmpar_exch_halo_field(field_4d_real)
+        ! Set up OpenMP threading.
+        call self % debug_print('Setting up OpenMP threading')
 
-                        nullify(field_4d_real)
-                    case (5)
-                        call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, &
-                            trim(adjustl(field_name)), field_5d_real, timelevel=1)
+        call mpas_atm_threading_init(self % domain_ptr % blocklist, ierr=ierr)
 
-                        if (.not. associated(field_5d_real)) then
-                            call self % model_error('Failed to find field "' // trim(adjustl(field_name)) // &
-                                '"', subname, __LINE__)
-                        end if
+        if (ierr /= 0) then
+            call self % model_error('OpenMP threading setup failed for core ' // trim(self % domain_ptr % core % corename), &
+                subname, __LINE__)
+        end if
 
-                        call mpas_dmpar_exch_halo_field(field_5d_real)
+        ! Set up inner dimensions used by arrays in optimized dynamics subroutines.
+        call self % debug_print('Setting up dimensions')
 
-                        nullify(field_5d_real)
-                    case default
-                        call self % model_error('Unsupported field rank ' // stringify([mpas_pool_field_info % ndims]), &
-                            subname, __LINE__)
-                end select
-            case default
-                call self % model_error('Unsupported field type (Must be one of: integer, real)', subname, __LINE__)
-        end select
+        call self % get_variable_pointer(nvertlevels, 'dim', 'nVertLevels')
+        call self % get_variable_pointer(maxedges, 'dim', 'maxEdges')
+        call self % get_variable_pointer(maxedges2, 'dim', 'maxEdges2')
+        call self % get_variable_pointer(num_scalars, 'dim', 'num_scalars')
+
+        call mpas_atm_set_dims(nvertlevels, maxedges, maxedges2, num_scalars)
+
+        nullify(nvertlevels, maxedges, maxedges2, num_scalars)
+
+        ! Build halo exchange groups and set the `exchange_halo_group` procedure pointer, which is used to
+        ! exchange the halo layers of all fields in the named group.
+        nullify(exchange_halo_group)
+
+        call atm_build_halo_groups(self % domain_ptr, ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to build halo exchange groups', subname, __LINE__)
+        end if
+
+        if (.not. associated(exchange_halo_group)) then
+            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
+        ! such variables into all subsequent time levels to initialize them.
+        call self % get_variable_pointer(config_do_restart, 'cfg', 'config_do_restart')
+
+        if (.not. config_do_restart) then
+            ! Run type is initial run.
+            call self % debug_print('Initializing time levels')
+
+            call self % get_pool_pointer(mpas_pool, 'state')
+
+            call mpas_pool_initialize_time_levels(mpas_pool)
+
+            nullify(mpas_pool)
+        end if
+
+        nullify(config_do_restart)
+
+        call exchange_halo_group(self % domain_ptr, 'initialization:u', ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to exchange halo layers for group "initialization:u"', subname, __LINE__)
+        end if
+
+        ! Initialize atmospheric variables (e.g., momentum, thermodynamic... variables in governing equations)
+        ! as well as various aspects of time in MPAS.
+
+        call self % debug_print('Initializing atmospheric variables')
+
+        ! Controlled by `config_start_time` in namelist.
+        mpas_time = mpas_get_clock_time(self % domain_ptr % clock, mpas_start_time, ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__)
+        end if
+
+        call mpas_get_time(mpas_time, datetimestring=date_time, ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__)
+        end if
+
+        ! Controlled by `config_dt` in namelist.
+        call self % get_pool_pointer(mpas_pool, 'mesh')
+        call self % get_variable_pointer(config_dt, 'cfg', 'config_dt')
+
+        call atm_mpas_init_block(self % domain_ptr % dminfo, self % domain_ptr % streammanager, self % domain_ptr % blocklist, &
+            mpas_pool, config_dt)
+
+        nullify(mpas_pool)
+        nullify(config_dt)
+
+        call self % get_variable_pointer(xtime, 'state', 'xtime', time_level=1)
+
+        xtime = date_time
+
+        nullify(xtime)
+
+        ! Initialize `initial_time` in the second time level. We need to do this manually because initial states
+        ! are read into time level 1, and if we write anything from time level 2, `initial_time` will be invalid.
+        call self % get_variable_pointer(initial_time_1, 'state', 'initial_time', time_level=1)
+        call self % get_variable_pointer(initial_time_2, 'state', 'initial_time', time_level=2)
+
+        initial_time_2 = initial_time_1
+
+        ! Set time units to CF-compliant "seconds since <date and time>".
+        call self % get_pool_pointer(mpas_pool, 'state')
+
+        call mpas_pool_get_field(mpas_pool, 'Time', field_0d_real, timelevel=1)
+
+        if (.not. associated(field_0d_real)) then
+            call self % model_error('Failed to find variable "Time"', subname, __LINE__)
+        end if
+
+        call mpas_modify_att(field_0d_real % attlists(1) % attlist, 'units', &
+            'seconds since ' // mpas_string_replace(initial_time_1, '_', ' '), ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to set time units', subname, __LINE__)
+        end if
+
+        nullify(initial_time_1, initial_time_2)
+        nullify(mpas_pool)
+        nullify(field_0d_real)
+
+        call exchange_halo_group(self % domain_ptr, 'initialization:pv_edge,ru,rw', ierr=ierr)
+
+        if (ierr /= 0) then
+            call self % model_error('Failed to exchange halo layers for group "initialization:pv_edge,ru,rw"', subname, __LINE__)
+        end if
+
+        call self % debug_print('Initializing dynamics')
+
+        ! Prepare dynamics for time integration.
+        call mpas_atm_dynamics_init(self % domain_ptr)
 
         call self % debug_print(subname // ' completed')
-    end subroutine dyn_mpas_exchange_halo
+
+        call self % debug_print('Successful initialization of MPAS dynamical core')
+    contains
+        !> Test if `a` is divisible by `b`, where `a` and `b` are both reals.
+        !> (KCW, 2024-05-25)
+        pure function almost_divisible(a, b)
+            real(rkind), intent(in) :: a, b
+            logical :: almost_divisible
+
+            real(rkind) :: error_tolerance
+
+            error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b))
+
+            if (almost_equal(mod(abs(a), abs(b)), 0.0_rkind, absolute_tolerance=error_tolerance) .or. &
+                almost_equal(mod(abs(a), abs(b)), abs(b), absolute_tolerance=error_tolerance)) then
+                almost_divisible = .true.
+
+                return
+            end if
+
+            almost_divisible = .false.
+        end function almost_divisible
+
+        !> Test `a` and `b` for approximate equality, where `a` and `b` are both reals.
+        !> (KCW, 2024-05-25)
+        pure function almost_equal(a, b, absolute_tolerance, relative_tolerance)
+            real(rkind), intent(in) :: a, b
+            real(rkind), optional, intent(in) :: absolute_tolerance, relative_tolerance
+            logical :: almost_equal
+
+            real(rkind) :: error_tolerance
+
+            if (present(relative_tolerance)) then
+                error_tolerance = relative_tolerance * max(abs(a), abs(b))
+            else
+                error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b))
+            end if
+
+            if (present(absolute_tolerance)) then
+                error_tolerance = max(absolute_tolerance, error_tolerance)
+            end if
+
+            if (abs(a - b) <= error_tolerance) then
+                almost_equal = .true.
+
+                return
+            end if
+
+            almost_equal = .false.
+        end function almost_equal
+    end subroutine dyn_mpas_init_phase4
+
+    !-------------------------------------------------------------------------------
+    ! function dyn_mpas_get_constituent_name
+    !
+    !> \brief  Query constituent name by its index
+    !> \author Kuan-Chih Wang
+    !> \date   2024-05-16
+    !> \details
+    !>  This function returns the constituent name that corresponds to the given
+    !>  constituent index. In case of errors, an empty character string is produced.
+    !
+    !-------------------------------------------------------------------------------
+    pure function dyn_mpas_get_constituent_name(self, constituent_index) result(constituent_name)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        integer, intent(in) :: constituent_index
+
+        character(:), allocatable :: constituent_name
+
+        ! Catch segmentation fault.
+        if (.not. allocated(self % constituent_name)) then
+            constituent_name = ''
+
+            return
+        end if
+
+        if (constituent_index < lbound(self % constituent_name, 1) .or. &
+            constituent_index > ubound(self % constituent_name, 1)) then
+            constituent_name = ''
+
+            return
+        end if
+
+        constituent_name = trim(adjustl(self % constituent_name(constituent_index)))
+    end function dyn_mpas_get_constituent_name
+
+    !-------------------------------------------------------------------------------
+    ! function dyn_mpas_get_constituent_index
+    !
+    !> \brief  Query constituent index by its name
+    !> \author Kuan-Chih Wang
+    !> \date   2024-05-16
+    !> \details
+    !>  This function returns the constituent index that corresponds to the given
+    !>  constituent name. In case of errors, zero is produced.
+    !
+    !-------------------------------------------------------------------------------
+    pure function dyn_mpas_get_constituent_index(self, constituent_name) result(constituent_index)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        character(*), intent(in) :: constituent_name
+
+        integer :: i
+        integer :: constituent_index
+
+        ! Catch segmentation fault.
+        if (.not. allocated(self % constituent_name)) then
+            constituent_index = 0
+
+            return
+        end if
+
+        do i = 1, self % number_of_constituents
+            if (trim(adjustl(constituent_name)) == trim(adjustl(self % constituent_name(i)))) then
+                constituent_index = i
+
+                return
+            end if
+        end do
+
+        constituent_index = 0
+    end function dyn_mpas_get_constituent_index
+
+    !-------------------------------------------------------------------------------
+    ! function dyn_mpas_map_mpas_scalar_index
+    !
+    !> \brief  Map MPAS scalar index from constituent index
+    !> \author Kuan-Chih Wang
+    !> \date   2024-05-16
+    !> \details
+    !>  This function returns the MPAS scalar index that corresponds to the given
+    !>  constituent index. In case of errors, zero is produced.
+    !
+    !-------------------------------------------------------------------------------
+    pure function dyn_mpas_map_mpas_scalar_index(self, constituent_index) result(mpas_scalar_index)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        integer, intent(in) :: constituent_index
+
+        integer :: mpas_scalar_index
+
+        ! Catch segmentation fault.
+        if (.not. allocated(self % index_constituent_to_mpas_scalar)) then
+            mpas_scalar_index = 0
+
+            return
+        end if
+
+        if (constituent_index < lbound(self % index_constituent_to_mpas_scalar, 1) .or. &
+            constituent_index > ubound(self % index_constituent_to_mpas_scalar, 1)) then
+            mpas_scalar_index = 0
+
+            return
+        end if
+
+        mpas_scalar_index = self % index_constituent_to_mpas_scalar(constituent_index)
+    end function dyn_mpas_map_mpas_scalar_index
+
+    !-------------------------------------------------------------------------------
+    ! function dyn_mpas_map_constituent_index
+    !
+    !> \brief  Map constituent index from MPAS scalar index
+    !> \author Kuan-Chih Wang
+    !> \date   2024-05-16
+    !> \details
+    !>  This function returns the constituent index that corresponds to the given
+    !>  MPAS scalar index. In case of errors, zero is produced.
+    !
+    !-------------------------------------------------------------------------------
+    pure function dyn_mpas_map_constituent_index(self, mpas_scalar_index) result(constituent_index)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        integer, intent(in) :: mpas_scalar_index
+
+        integer :: constituent_index
+
+        ! Catch segmentation fault.
+        if (.not. allocated(self % index_mpas_scalar_to_constituent)) then
+            constituent_index = 0
+
+            return
+        end if
+
+        if (mpas_scalar_index < lbound(self % index_mpas_scalar_to_constituent, 1) .or. &
+            mpas_scalar_index > ubound(self % index_mpas_scalar_to_constituent, 1)) then
+            constituent_index = 0
+
+            return
+        end if
+
+        constituent_index = self % index_mpas_scalar_to_constituent(mpas_scalar_index)
+    end function dyn_mpas_map_constituent_index
+
+    !-------------------------------------------------------------------------------
+    ! subroutine dyn_mpas_get_local_mesh_dimension
+    !
+    !> \brief  Returns 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
+    !>    on each individual task, both with/without halo points.
+    !
+    !-------------------------------------------------------------------------------
+    subroutine dyn_mpas_get_local_mesh_dimension(self, &
+            ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels)
+        class(mpas_dynamical_core_type), intent(in) :: self
+        integer, intent(out) :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
+
+        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_local_mesh_dimension'
+        integer, pointer :: ncells_pointer
+        integer, pointer :: ncellssolve_pointer
+        integer, pointer :: nedges_pointer
+        integer, pointer :: nedgessolve_pointer
+        integer, pointer :: nvertices_pointer
+        integer, pointer :: nverticessolve_pointer
+        integer, pointer :: nvertlevels_pointer
+
+        nullify(ncells_pointer)
+        nullify(ncellssolve_pointer)
+        nullify(nedges_pointer)
+        nullify(nedgessolve_pointer)
+        nullify(nvertices_pointer)
+        nullify(nverticessolve_pointer)
+        nullify(nvertlevels_pointer)
+
+        call self % get_variable_pointer(ncells_pointer, 'dim', 'nCells')
+        call self % get_variable_pointer(ncellssolve_pointer, 'dim', 'nCellsSolve')
+        call self % get_variable_pointer(nedges_pointer, 'dim', 'nEdges')
+        call self % get_variable_pointer(nedgessolve_pointer, 'dim', 'nEdgesSolve')
+        call self % get_variable_pointer(nvertices_pointer, 'dim', 'nVertices')
+        call self % get_variable_pointer(nverticessolve_pointer, 'dim', 'nVerticesSolve')
+        call self % get_variable_pointer(nvertlevels_pointer, 'dim', 'nVertLevels')
+
+        ncells = ncells_pointer                  ! Number of cells, including halo cells.
+        ncells_solve = ncellssolve_pointer       ! Number of cells, excluding halo cells.
+        nedges = nedges_pointer                  ! Number of edges, including halo edges.
+        nedges_solve = nedgessolve_pointer       ! Number of edges, excluding halo edges.
+        nvertices = nvertices_pointer            ! Number of vertices, including halo vertices.
+        nvertices_solve = nverticessolve_pointer ! Number of vertices, excluding halo vertices.
+
+        ! Vertical levels are not decomposed.
+        ! All tasks have the same number of vertical levels.
+        nvertlevels = nvertlevels_pointer
+
+        nullify(ncells_pointer)
+        nullify(ncellssolve_pointer)
+        nullify(nedges_pointer)
+        nullify(nedgessolve_pointer)
+        nullify(nvertices_pointer)
+        nullify(nverticessolve_pointer)
+        nullify(nvertlevels_pointer)
+    end subroutine dyn_mpas_get_local_mesh_dimension
 
     !-------------------------------------------------------------------------------
     ! subroutine dyn_mpas_get_global_mesh_dimension
@@ -1597,11 +3067,17 @@ subroutine dyn_mpas_get_global_mesh_dimension(self, &
         real(rkind), intent(out) :: sphere_radius
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_global_mesh_dimension'
-        integer, pointer :: maxedges_pointer => null()
-        integer, pointer :: ncellssolve_pointer => null()
-        integer, pointer :: nedgessolve_pointer => null()
-        integer, pointer :: nverticessolve_pointer => null()
-        integer, pointer :: nvertlevels_pointer => null()
+        integer, pointer :: maxedges_pointer
+        integer, pointer :: ncellssolve_pointer
+        integer, pointer :: nedgessolve_pointer
+        integer, pointer :: nverticessolve_pointer
+        integer, pointer :: nvertlevels_pointer
+
+        nullify(maxedges_pointer)
+        nullify(ncellssolve_pointer)
+        nullify(nedgessolve_pointer)
+        nullify(nverticessolve_pointer)
+        nullify(nvertlevels_pointer)
 
         call self % get_variable_pointer(maxedges_pointer, 'dim', 'maxEdges')
         call self % get_variable_pointer(ncellssolve_pointer, 'dim', 'nCellsSolve')
@@ -1685,8 +3161,9 @@ subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_c0'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
 
@@ -1712,8 +3189,9 @@ subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_c1'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1733,8 +3211,9 @@ subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i0'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
 
@@ -1763,8 +3242,9 @@ subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i1'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
 
@@ -1790,8 +3270,9 @@ subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i2'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1811,8 +3292,9 @@ subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i3'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1832,8 +3314,9 @@ subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_l0'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
 
@@ -1857,8 +3340,9 @@ subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r0'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
 
@@ -1884,8 +3368,9 @@ subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r1'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1905,8 +3390,9 @@ subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r2'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1926,8 +3412,9 @@ subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r3'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1947,8 +3434,9 @@ subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r4'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -1968,8 +3456,9 @@ subroutine dyn_mpas_get_variable_pointer_r5(self, variable_pointer, pool_name, v
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r5'
-        type(mpas_pool_type), pointer :: mpas_pool => null()
+        type(mpas_pool_type), pointer :: mpas_pool
 
+        nullify(mpas_pool)
         call self % get_pool_pointer(mpas_pool, pool_name)
         nullify(variable_pointer)
         call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level)
@@ -2002,9 +3491,10 @@ subroutine dyn_mpas_get_variable_value_c0(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_c0'
-        character(strkind), pointer :: variable_pointer => null()
+        character(strkind), pointer :: variable_pointer
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2023,9 +3513,10 @@ subroutine dyn_mpas_get_variable_value_c1(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_c1'
-        character(strkind), pointer :: variable_pointer(:) => null()
+        character(strkind), pointer :: variable_pointer(:)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2044,9 +3535,10 @@ subroutine dyn_mpas_get_variable_value_i0(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i0'
-        integer, pointer :: variable_pointer => null()
+        integer, pointer :: variable_pointer
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2065,9 +3557,10 @@ subroutine dyn_mpas_get_variable_value_i1(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i1'
-        integer, pointer :: variable_pointer(:) => null()
+        integer, pointer :: variable_pointer(:)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2086,9 +3579,10 @@ subroutine dyn_mpas_get_variable_value_i2(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i2'
-        integer, pointer :: variable_pointer(:, :) => null()
+        integer, pointer :: variable_pointer(:, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2107,9 +3601,10 @@ subroutine dyn_mpas_get_variable_value_i3(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i3'
-        integer, pointer :: variable_pointer(:, :, :) => null()
+        integer, pointer :: variable_pointer(:, :, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2128,9 +3623,10 @@ subroutine dyn_mpas_get_variable_value_l0(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_l0'
-        logical, pointer :: variable_pointer => null()
+        logical, pointer :: variable_pointer
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2149,9 +3645,10 @@ subroutine dyn_mpas_get_variable_value_r0(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r0'
-        real(rkind), pointer :: variable_pointer => null()
+        real(rkind), pointer :: variable_pointer
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2170,9 +3667,10 @@ subroutine dyn_mpas_get_variable_value_r1(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r1'
-        real(rkind), pointer :: variable_pointer(:) => null()
+        real(rkind), pointer :: variable_pointer(:)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2191,9 +3689,10 @@ subroutine dyn_mpas_get_variable_value_r2(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r2'
-        real(rkind), pointer :: variable_pointer(:, :) => null()
+        real(rkind), pointer :: variable_pointer(:, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2212,9 +3711,10 @@ subroutine dyn_mpas_get_variable_value_r3(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r3'
-        real(rkind), pointer :: variable_pointer(:, :, :) => null()
+        real(rkind), pointer :: variable_pointer(:, :, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2233,9 +3733,10 @@ subroutine dyn_mpas_get_variable_value_r4(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r4'
-        real(rkind), pointer :: variable_pointer(:, :, :, :) => null()
+        real(rkind), pointer :: variable_pointer(:, :, :, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
@@ -2254,9 +3755,10 @@ subroutine dyn_mpas_get_variable_value_r5(self, variable_value, pool_name, varia
         integer, optional, intent(in) :: time_level
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r5'
-        real(rkind), pointer :: variable_pointer(:, :, :, :, :) => null()
+        real(rkind), pointer :: variable_pointer(:, :, :, :, :)
         integer :: ierr
 
+        nullify(variable_pointer)
         call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level)
         allocate(variable_value, source=variable_pointer, stat=ierr)
 
diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90
index 5b4af9cd..5c2d8f0b 100644
--- a/src/dynamics/mpas/dyn_comp.F90
+++ b/src/dynamics/mpas/dyn_comp.F90
@@ -1,28 +1,50 @@
 module dyn_comp
     use dyn_mpas_subdriver, only: mpas_dynamical_core_type
 
-    ! Modules from CAM.
-    use cam_abortutils, only: endrun
+    ! Modules from CAM-SIMA.
+    use air_composition, only: thermodynamic_active_species_num, &
+                               thermodynamic_active_species_liq_num, &
+                               thermodynamic_active_species_ice_num, &
+                               thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore, &
+                               thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, &
+                               thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore
+    use cam_abortutils, only: check_allocate, endrun
+    use cam_constituents, only: const_name, const_is_water_species, num_advected, readtrace
     use cam_control_mod, only: initial_run
+    use cam_field_read, only: cam_read_field
+    use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id
+    use cam_initfiles, only: initial_file_get_id, topo_file_get_id
     use cam_instance, only: atm_id
     use cam_logfile, only: debug_output, debugout_none, iulog
+    use cam_pio_utils, only: clean_iodesc_list
+    use dyn_tests_utils, only: vc_height
+    use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, constant_pi => pi, &
+                        constant_rd => rair, constant_rv => rh2o, &
+                        deg_to_rad
+    use inic_analytic, only: analytic_ic_active, dyn_set_inic_col
     use runtime_obj, only: runtime_options
     use spmd_utils, only: iam, masterproc, mpicom
     use string_utils, only: stringify
-    use time_manager, only: get_start_date, get_stop_date, get_run_duration, timemgr_get_calendar_cf
+    use time_manager, only: get_start_date, get_stop_date, get_step_size, get_run_duration, timemgr_get_calendar_cf
+    use vert_coord, only: pver, pverp
 
     ! Modules from CESM Share.
     use shr_file_mod, only: shr_file_getunit
-    use shr_kind_mod, only: shr_kind_cs
+    use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8
     use shr_pio_mod, only: shr_pio_getiosys
 
+    ! Modules from CCPP.
+    use cam_ccpp_cap, only: cam_constituents_array
+    use ccpp_kinds, only: kind_phys
+    use phys_vars_init_check, only: mark_as_initialized, std_name_len
+
     ! Modules from external libraries.
-    use pio, only: iosystem_desc_t
+    use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open
 
     implicit none
 
     private
-    ! Provide APIs required by CAM Control.
+    ! Provide APIs required by CAM-SIMA.
     public :: dyn_import_t
     public :: dyn_export_t
     public :: dyn_readnl
@@ -32,15 +54,31 @@ module dyn_comp
 
     public :: dyn_debug_print
     public :: mpas_dynamical_core
+    public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
+    public :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max
+    public :: sphere_radius
 
+    !> NOTE:
+    !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM-SIMA requires it.
+    !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object"
+    !> below.
     type :: dyn_import_t
     end type dyn_import_t
 
+    !> NOTE:
+    !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM-SIMA requires it.
+    !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object"
+    !> below.
     type :: dyn_export_t
     end type dyn_export_t
 
     !> The "instance/object" of MPAS dynamical core.
     type(mpas_dynamical_core_type) :: mpas_dynamical_core
+
+    ! Local and global mesh dimensions of MPAS dynamical core.
+    integer :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
+    integer :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max
+    real(kind_r8) :: sphere_radius
 contains
     !> Print a debug message with optionally the value(s) of a variable.
     !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead.
@@ -83,13 +121,15 @@ subroutine dyn_readnl(namelist_path)
         character(*), intent(in) :: namelist_path
 
         character(*), parameter :: subname = 'dyn_comp::dyn_readnl'
-        character(shr_kind_cs) :: cam_calendar
+        character(kind_cs) :: cam_calendar
         integer :: log_unit(2)
         integer :: start_date_time(6), & ! YYYY, MM, DD, hh, mm, ss.
                    stop_date_time(6),  & ! YYYY, MM, DD, hh, mm, ss.
                    run_duration(4),    & ! DD, hh, mm, ss.
                    sec_since_midnight    ! Second(s) since midnight.
-        type(iosystem_desc_t), pointer :: pio_iosystem => null()
+        type(iosystem_desc_t), pointer :: pio_iosystem
+
+        nullify(pio_iosystem)
 
         ! Enable/disable the debug output of MPAS dynamical core according to the debug verbosity level of CAM-SIMA.
         mpas_dynamical_core % debug_output = (debug_output > debugout_none)
@@ -99,7 +139,6 @@ subroutine dyn_readnl(namelist_path)
         log_unit(2) = shr_file_getunit()
 
         ! Initialize MPAS framework with supplied MPI communicator group and log units.
-        ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details.
         call mpas_dynamical_core % init_phase1(mpicom, endrun, iulog, log_unit)
 
         cam_calendar = timemgr_get_calendar_cf()
@@ -114,14 +153,12 @@ subroutine dyn_readnl(namelist_path)
         run_duration(2:4) = sec_to_hour_min_sec(sec_since_midnight)
 
         ! Read MPAS-related namelist variables from `namelist_path`, e.g., `atm_in`.
-        ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details.
         call mpas_dynamical_core % read_namelist(namelist_path, &
             cam_calendar, start_date_time, stop_date_time, run_duration, initial_run)
 
         pio_iosystem => shr_pio_getiosys(atm_id)
 
         ! Initialize MPAS framework with supplied PIO system descriptor.
-        ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details.
         call mpas_dynamical_core % init_phase2(pio_iosystem)
 
         nullify(pio_iosystem)
@@ -139,13 +176,664 @@ pure function sec_to_hour_min_sec(sec) result(hour_min_sec)
         end function sec_to_hour_min_sec
     end subroutine dyn_readnl
 
+    !> Initialize MPAS dynamical core by one of the following:
+    !> 1. Setting analytic initial condition;
+    !> 2. Reading initial condition from a file;
+    !> 3. Restarting from a file.
+    !> (KCW, 2024-05-28)
+    !
     ! Called by `cam_init` in `src/control/cam_comp.F90`.
     subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
-        type(runtime_options), intent(in)  :: cam_runtime_opts
-        type(dyn_import_t),    intent(out) :: dyn_in
-        type(dyn_export_t),    intent(out) :: dyn_out
+        type(runtime_options), intent(in) :: cam_runtime_opts
+        type(dyn_import_t), intent(in) :: dyn_in
+        type(dyn_export_t), intent(in) :: dyn_out
+
+        character(*), parameter :: subname = 'dyn_comp::dyn_init'
+        character(std_name_len), allocatable :: constituent_name(:)
+        integer :: coupling_time_interval
+        integer :: i
+        integer :: ierr
+        logical, allocatable :: is_water_species(:)
+        type(file_desc_t), pointer :: pio_init_file
+        type(file_desc_t), pointer :: pio_topo_file
+
+        nullify(pio_init_file)
+        nullify(pio_topo_file)
+
+        allocate(constituent_name(num_advected), stat=ierr)
+        call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__)
+
+        allocate(is_water_species(num_advected), stat=ierr)
+        call check_allocate(ierr, subname, 'is_water_species(num_advected)', 'dyn_comp', __LINE__)
+
+        do i = 1, num_advected
+            constituent_name(i) = const_name(i)
+            is_water_species(i) = const_is_water_species(i)
+        end do
+
+        ! Inform MPAS about constituent names and their corresponding waterness.
+        call mpas_dynamical_core % define_scalar(constituent_name, is_water_species)
+
+        deallocate(constituent_name)
+        deallocate(is_water_species)
+
+        ! Provide mapping information between MPAS scalars and constituent names to CAM-SIMA.
+        do i = 1, thermodynamic_active_species_num
+            thermodynamic_active_species_idx_dycore(i) = &
+                mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_idx(i))
+        end do
+
+        do i = 1, thermodynamic_active_species_liq_num
+            thermodynamic_active_species_liq_idx_dycore(i) = &
+                mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_liq_idx(i))
+        end do
+
+        do i = 1, thermodynamic_active_species_ice_num
+            thermodynamic_active_species_ice_idx_dycore(i) = &
+                mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_ice_idx(i))
+        end do
+
+        pio_init_file => initial_file_get_id()
+        pio_topo_file => topo_file_get_id()
+
+        if (initial_run) then
+            ! Run type is initial run.
+
+            call dyn_debug_print('Calling check_topography_data')
+
+            call check_topography_data(pio_topo_file)
+
+            if (analytic_ic_active()) then
+                call dyn_debug_print('Calling set_analytic_initial_condition')
+
+                call set_analytic_initial_condition()
+            else
+                ! Perform default initialization for all constituents.
+                ! Subsequently, they can be overridden depending on the namelist option (below) and
+                ! the actual availability (checked and handled by MPAS).
+                call dyn_debug_print('Calling set_default_constituent')
+
+                call set_default_constituent()
+
+                ! Namelist option that controls if constituents are to be read from the file.
+                if (readtrace) then
+                    ! Read variables that belong to the "input" stream in MPAS.
+                    call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input')
+                else
+                    ! Read variables that belong to the "input" stream in MPAS, excluding constituents.
+                    call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input-scalars')
+                end if
+            end if
+        else
+            ! Run type is branch or restart run.
+
+            ! Read variables that belong to the "input" and "restart" streams in MPAS.
+            call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input+restart')
+        end if
+
+        call clean_iodesc_list()
+        call mark_variable_as_initialized()
+
+        nullify(pio_init_file)
+        nullify(pio_topo_file)
+
+        ! This is the time interval for dynamics-physics coupling in CAM-SIMA.
+        ! Each time MPAS dynamical core is called to run, it will integrate with time for this specific interval,
+        ! then yield control back to the caller.
+        coupling_time_interval = get_step_size()
+
+        ! Finish MPAS dynamical core initialization. After this point, MPAS dynamical core is ready for time integration.
+        call mpas_dynamical_core % init_phase4(coupling_time_interval)
     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,
+    !> 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)
+    subroutine check_topography_data(pio_file)
+        type(file_desc_t), pointer, intent(in) :: pio_file
+
+        character(*), parameter :: subname = 'dyn_comp::check_topography_data'
+        integer :: ierr
+        logical :: success
+        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_r8), pointer :: zgrid(:, :)                        ! From MPAS. Geometric height (meters) at layer interfaces.
+
+        nullify(zgrid)
+
+        call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
+
+        if (associated(pio_file)) then
+            call dyn_debug_print('Topography file is used')
+
+            if (.not. pio_file_is_open(pio_file)) then
+                call endrun('Invalid PIO file descriptor', subname, __LINE__)
+            end if
+
+            allocate(surface_geopotential(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'surface_geopotential(ncells_solve)', 'dyn_comp', __LINE__)
+
+            allocate(surface_geometric_height(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'surface_geometric_height(ncells_solve)', 'dyn_comp', __LINE__)
+
+            surface_geopotential(:) = 0.0_kind_r8
+            surface_geometric_height(:) = 0.0_kind_r8
+
+            call cam_read_field('PHIS', pio_file, surface_geopotential, success, &
+                gridname='cam_cell', timelevel=1, log_output=(debug_output > debugout_none))
+
+            if (.not. success) then
+                call endrun('Failed to find variable "PHIS"', subname, __LINE__)
+            end if
+
+            surface_geometric_height(:) = surface_geopotential(:) / constant_g
+
+            ! Surface geometric height in MPAS should match the values in topography file.
+            if (any(abs(zgrid(1, 1:ncells_solve) - surface_geometric_height) > error_tolerance)) then
+                call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__)
+            end if
+
+            deallocate(surface_geopotential)
+            deallocate(surface_geometric_height)
+        else
+            call dyn_debug_print('Topography file is not used')
+
+            ! Surface geometric height in MPAS should be zero.
+            if (any(abs(zgrid(1, 1:ncells_solve)) > error_tolerance)) then
+                call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__)
+            end if
+        end if
+
+        nullify(zgrid)
+    end subroutine check_topography_data
+
+    !> Set analytic initial condition for MPAS.
+    !> (KCW, 2024-05-22)
+    subroutine set_analytic_initial_condition()
+        character(*), parameter :: subname = 'dyn_comp::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.
+                                                        ! Dimension and vertical index orders follow CAM-SIMA convention.
+        real(kind_r8), pointer :: zgrid(:, :)           ! Geometric height (meters) at layer interfaces.
+                                                        ! Dimension and vertical index orders follow MPAS convention.
+
+        call init_shared_variable()
+
+        call set_mpas_state_u()
+        call set_mpas_state_w()
+        call set_mpas_state_scalars()
+        call set_mpas_state_rho_theta()
+        call set_mpas_state_rho_base_theta_base()
+
+        deallocate(global_grid_index)
+        deallocate(lat_rad, lon_rad)
+        deallocate(z_int)
+
+        nullify(zgrid)
+    contains
+        !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines.
+        !> (KCW, 2024-05-13)
+        subroutine init_shared_variable()
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variable'
+            integer :: ierr
+            integer :: k
+            integer, pointer :: indextocellid(:)
+            real(kind_r8), pointer :: lat_deg(:), lon_deg(:)
+
+            call dyn_debug_print('Preparing to set analytic initial condition')
+
+            nullify(zgrid)
+            nullify(indextocellid)
+            nullify(lat_deg, lon_deg)
+
+            allocate(global_grid_index(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'global_grid_index(ncells_solve)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(indextocellid, 'mesh', 'indexToCellID')
+
+            global_grid_index(:) = indextocellid(1:ncells_solve)
+
+            nullify(indextocellid)
+
+            allocate(lat_rad(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'lat_rad(ncells_solve)', 'dyn_comp', __LINE__)
+
+            allocate(lon_rad(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'lon_rad(ncells_solve)', 'dyn_comp', __LINE__)
+
+            ! "mpas_cell" is a registered grid name that is defined in `dyn_grid`.
+            lat_deg => cam_grid_get_latvals(cam_grid_id('mpas_cell'))
+            lon_deg => cam_grid_get_lonvals(cam_grid_id('mpas_cell'))
+
+            if (.not. associated(lat_deg)) then
+                call endrun('Failed to find variable "lat_deg"', subname, __LINE__)
+            end if
+
+            if (.not. associated(lon_deg)) then
+                call endrun('Failed to find variable "lon_deg"', subname, __LINE__)
+            end if
+
+            lat_rad(:) = lat_deg(:) * deg_to_rad
+            lon_rad(:) = lon_deg(:) * deg_to_rad
+
+            nullify(lat_deg, lon_deg)
+
+            allocate(z_int(ncells_solve, pverp), stat=ierr)
+            call check_allocate(ierr, subname, 'z_int(ncells_solve, pverp)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
+
+            ! Vertical index order is reversed between CAM-SIMA and MPAS.
+            do k = 1, pverp
+                z_int(:, k) = zgrid(pverp - k + 1, 1:ncells_solve)
+            end do
+        end subroutine init_shared_variable
+
+        !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces).
+        !> (KCW, 2024-05-13)
+        subroutine set_mpas_state_u()
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u'
+            integer :: ierr
+            integer :: k
+            real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :)
+
+            call dyn_debug_print('Setting MPAS state "u"')
+
+            nullify(ucellzonal, ucellmeridional)
+
+            allocate(buffer_2d_real(ncells_solve, pver), stat=ierr)
+            call check_allocate(ierr, subname, 'buffer_2d_real(ncells_solve, pver)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
+            call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
+
+            buffer_2d_real(:, :) = 0.0_kind_r8
+
+            call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, u=buffer_2d_real)
+
+            ! Vertical index order is reversed between CAM-SIMA and MPAS.
+            do k = 1, pver
+                ucellzonal(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1)
+            end do
+
+            buffer_2d_real(:, :) = 0.0_kind_r8
+
+            call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, v=buffer_2d_real)
+
+            ! Vertical index order is reversed between CAM-SIMA and MPAS.
+            do k = 1, pver
+                ucellmeridional(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1)
+            end do
+
+            deallocate(buffer_2d_real)
+
+            nullify(ucellzonal, ucellmeridional)
+
+            call mpas_dynamical_core % compute_edge_wind()
+        end subroutine set_mpas_state_u
+
+        !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces).
+        !> (KCW, 2024-05-13)
+        subroutine set_mpas_state_w()
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_w'
+            real(kind_r8), pointer :: w(:, :)
+
+            call dyn_debug_print('Setting MPAS state "w"')
+
+            nullify(w)
+
+            call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1)
+
+            w(:, 1:ncells_solve) = 0.0_kind_r8
+
+            nullify(w)
+
+            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+            call mpas_dynamical_core % exchange_halo('w')
+        end subroutine set_mpas_state_w
+
+        !> Set MPAS state `scalars` (i.e., constituents).
+        !> (KCW, 2024-05-17)
+        subroutine set_mpas_state_scalars()
+            ! CCPP standard name of `qv`, which denotes water vapor mixing ratio.
+            character(*), parameter :: constituent_qv_standard_name = &
+                'water_vapor_mixing_ratio_wrt_dry_air'
+
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_scalars'
+            integer :: i, k
+            integer :: ierr
+            integer, allocatable :: constituent_index(:)
+            integer, pointer :: index_qv
+            real(kind_r8), pointer :: scalars(:, :, :)
+
+            call dyn_debug_print('Setting MPAS state "scalars"')
+
+            nullify(index_qv)
+            nullify(scalars)
+
+            allocate(buffer_3d_real(ncells_solve, pver, num_advected), stat=ierr)
+            call check_allocate(ierr, subname, 'buffer_3d_real(ncells_solve, pver, num_advected)', 'dyn_comp', __LINE__)
+
+            allocate(constituent_index(num_advected), stat=ierr)
+            call check_allocate(ierr, subname, 'constituent_index(num_advected)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
+            call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
+
+            buffer_3d_real(:, :, :) = 0.0_kind_r8
+            constituent_index(:) = [(i, i = 1, num_advected)]
+
+            call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, q=buffer_3d_real, &
+                m_cnst=constituent_index)
+
+            ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index.
+            do i = 1, num_advected
+                ! Vertical index order is reversed between CAM-SIMA and MPAS.
+                do k = 1, pver
+                    scalars(i, k, 1:ncells_solve) = &
+                        buffer_3d_real(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i))
+                end do
+            end do
+
+            if (mpas_dynamical_core % get_constituent_name(mpas_dynamical_core % map_constituent_index(index_qv)) == &
+                constituent_qv_standard_name) then
+                ! The definition of `qv` matches exactly what MPAS wants. No conversion is needed.
+                call dyn_debug_print('No conversion is needed for water vapor mixing ratio')
+            else
+                ! The definition of `qv` actually represents specific humidity. Conversion is needed.
+                call dyn_debug_print('Conversion is needed and applied for water vapor mixing ratio')
+
+                ! Convert specific humidity to water vapor mixing ratio.
+                scalars(index_qv, :, 1:ncells_solve) = &
+                    scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_r8 - scalars(index_qv, :, 1:ncells_solve))
+            end if
+
+            deallocate(buffer_3d_real)
+            deallocate(constituent_index)
+
+            nullify(index_qv)
+            nullify(scalars)
+
+            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+            call mpas_dynamical_core % exchange_halo('scalars')
+        end subroutine set_mpas_state_scalars
+
+        !> Set MPAS state `rho` (i.e., dry air density) and `theta` (i.e., potential temperature).
+        !> (KCW, 2024-05-19)
+        subroutine set_mpas_state_rho_theta()
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_theta'
+            integer :: i, k
+            integer :: ierr
+            integer, pointer :: index_qv
+            real(kind_r8), allocatable :: p_mid_col(:)  ! Pressure (Pa) at layer midpoints of each column. This is full pressure,
+                                                        ! which also accounts for water vapor.
+            real(kind_r8), allocatable :: p_sfc(:)      ! Pressure (Pa) at surface. This is full pressure,
+                                                        ! which also accounts for water vapor.
+            real(kind_r8), allocatable :: qv_mid_col(:) ! Water vapor mixing ratio (kg/kg) at layer midpoints of each column.
+            real(kind_r8), allocatable :: t_mid(:, :)   ! Temperature (K) at layer midpoints.
+            real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column.
+            real(kind_r8), pointer :: rho(:, :)
+            real(kind_r8), pointer :: theta(:, :)
+            real(kind_r8), pointer :: scalars(:, :, :)
+
+            call dyn_debug_print('Setting MPAS state "rho" and "theta"')
+
+            nullify(index_qv)
+            nullify(rho)
+            nullify(theta)
+            nullify(scalars)
+
+            allocate(p_sfc(ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 'p_sfc(ncells_solve)', 'dyn_comp', __LINE__)
+
+            p_sfc(:) = 0.0_kind_r8
+
+            call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, ps=p_sfc)
+
+            allocate(buffer_2d_real(ncells_solve, pver), stat=ierr)
+            call check_allocate(ierr, subname, 'buffer_2d_real(ncells_solve, pver)', 'dyn_comp', __LINE__)
+
+            allocate(t_mid(pver, ncells_solve), stat=ierr)
+            call check_allocate(ierr, subname, 't_mid(pver, ncells_solve)', 'dyn_comp', __LINE__)
+
+            buffer_2d_real(:, :) = 0.0_kind_r8
+
+            call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, t=buffer_2d_real)
+
+            ! Vertical index order is reversed between CAM-SIMA and MPAS.
+            do k = 1, pver
+                t_mid(k, :) = buffer_2d_real(:, pver - k + 1)
+            end do
+
+            deallocate(buffer_2d_real)
+
+            allocate(p_mid_col(pver), stat=ierr)
+            call check_allocate(ierr, subname, 'p_mid_col(pver)', 'dyn_comp', __LINE__)
+
+            allocate(qv_mid_col(pver), stat=ierr)
+            call check_allocate(ierr, subname, 'qv_mid_col(pver)', 'dyn_comp', __LINE__)
+
+            allocate(tm_mid_col(pver), stat=ierr)
+            call check_allocate(ierr, subname, 'tm_mid_col(pver)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
+            call mpas_dynamical_core % get_variable_pointer(rho, 'diag', 'rho')
+            call mpas_dynamical_core % get_variable_pointer(theta, 'diag', 'theta')
+            call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
+
+            ! Set `rho` and `theta` column by column. This way, peak memory usage can be reduced.
+            do i = 1, ncells_solve
+                qv_mid_col(:) = scalars(index_qv, :, i)
+                tm_mid_col(:) = t_mid(:, i) * (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:))
+
+                ! Piecewise integrate hypsometric equation to derive `p_mid_col(1)`.
+                ! The formulation used here is exact.
+                p_mid_col(1) = p_by_hypsometric_equation( &
+                    p_sfc(i), &
+                    zgrid(1, i), &
+                    tm_mid_col(1) / (1.0_kind_r8 + qv_mid_col(1)), &
+                    0.5_kind_r8 * (zgrid(2, i) + zgrid(1, i)))
+
+                ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`.
+                ! The formulation used here is exact.
+                do k = 2, pver
+                    p_mid_col(k) = p_by_hypsometric_equation( &
+                        p_by_hypsometric_equation( &
+                            p_mid_col(k - 1), &
+                            0.5_kind_r8 * (zgrid(k, i) + zgrid(k - 1, i)), &
+                            tm_mid_col(k - 1) / (1.0_kind_r8 + qv_mid_col(k - 1)), &
+                            zgrid(k, i)), &
+                        zgrid(k, i), &
+                        tm_mid_col(k) / (1.0_kind_r8 + qv_mid_col(k)), &
+                        0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i)))
+                end do
+
+                rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:))
+                theta(:, i) = theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0)
+            end do
+
+            deallocate(p_mid_col)
+            deallocate(p_sfc)
+            deallocate(qv_mid_col)
+            deallocate(t_mid)
+            deallocate(tm_mid_col)
+
+            nullify(index_qv)
+            nullify(rho)
+            nullify(theta)
+            nullify(scalars)
+
+            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+            call mpas_dynamical_core % exchange_halo('rho')
+            call mpas_dynamical_core % exchange_halo('theta')
+        end subroutine set_mpas_state_rho_theta
+
+        !> Set MPAS state `rho_base` (i.e., base state dry air density) and `theta_base` (i.e., base state potential temperature).
+        !> (KCW, 2024-05-21)
+        subroutine set_mpas_state_rho_base_theta_base()
+            character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_base_theta_base'
+            integer :: i, k
+            integer :: ierr
+            real(kind_r8), parameter :: t_base = 250.0_kind_r8 ! Base state temperature (K) of dry isothermal atmosphere.
+                                                               ! The value used here is identical to MPAS.
+            real(kind_r8), allocatable :: p_base(:)            ! Base state pressure (Pa) at layer midpoints of each column.
+            real(kind_r8), pointer :: rho_base(:, :)
+            real(kind_r8), pointer :: theta_base(:, :)
+            real(kind_r8), pointer :: zz(:, :)
+
+            call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"')
+
+            nullify(rho_base)
+            nullify(theta_base)
+            nullify(zz)
+
+            allocate(p_base(pver), stat=ierr)
+            call check_allocate(ierr, subname, 'p_base(pver)', 'dyn_comp', __LINE__)
+
+            call mpas_dynamical_core % get_variable_pointer(rho_base, 'diag', 'rho_base')
+            call mpas_dynamical_core % get_variable_pointer(theta_base, 'diag', 'theta_base')
+            call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz')
+
+            ! Set `rho_base` and `theta_base` column by column. This way, peak memory usage can be reduced.
+            do i = 1, ncells_solve
+                do k = 1, pver
+                    ! Derive `p_base` by hypsometric equation.
+                    ! The formulation used here is exact and identical to MPAS.
+                    p_base(k) = p_by_hypsometric_equation( &
+                        constant_p0, &
+                        0.0_kind_r8, &
+                        t_base, &
+                        0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i)))
+                end do
+
+                rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i))
+                theta_base(:, i) = theta_by_poisson_equation(p_base, t_base, constant_p0)
+            end do
+
+            deallocate(p_base)
+
+            nullify(rho_base)
+            nullify(theta_base)
+            nullify(zz)
+
+            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+            call mpas_dynamical_core % exchange_halo('rho_base')
+            call mpas_dynamical_core % exchange_halo('theta_base')
+        end subroutine set_mpas_state_rho_base_theta_base
+
+        ! ----- p_2, z_2 ----- (Layer 2)
+        !       t_v
+        ! ----- p_1, z_1 ----- (Layer 1)
+        !
+        !> Compute the pressure `p_2` at height `z_2` from the pressure `p_1` at height `z_1` by hypsometric equation.
+        !> `t_v` is the mean virtual temperature between `z_1` and `z_2`.
+        !> (KCW, 2024-07-02)
+        pure elemental function p_by_hypsometric_equation(p_1, z_1, t_v, z_2) result(p_2)
+            real(kind_r8), intent(in) :: p_1, z_1, t_v, z_2
+            real(kind_r8) :: p_2
+
+            p_2 = p_1 * exp(-(z_2 - z_1) * constant_g / (constant_rd * t_v))
+        end function p_by_hypsometric_equation
+
+        ! ----- p_1, t_1 ----- (Arbitrary layer)
+        !
+        ! ----- p_0, t_0 ----- (Reference layer)
+        !
+        !> Compute the potential temperature `t_0` at reference pressure `p_0` from the temperature `t_1` at pressure `p_1` by
+        !> Poisson equation.
+        !> (KCW, 2024-07-02)
+        pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0)
+            real(kind_r8), intent(in) :: p_1, t_1, p_0
+            real(kind_r8) :: t_0
+
+            t_0 = t_1 * ((p_0 / p_1) ** (constant_rd / constant_cpd))
+        end function theta_by_poisson_equation
+    end subroutine set_analytic_initial_condition
+
+    !> Set default MPAS state `scalars` (i.e., constituents) in accordance with CCPP, which is a component of CAM-SIMA.
+    !> (KCW, 2024-07-09)
+    subroutine set_default_constituent()
+        character(*), parameter :: subname = 'dyn_comp::set_default_constituent'
+        integer :: i, k
+        real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory.
+        real(kind_r8), pointer :: scalars(:, :, :)        ! This points to MPAS memory.
+
+        call dyn_debug_print('Setting default MPAS state "scalars"')
+
+        nullify(constituents)
+        nullify(scalars)
+
+        constituents => cam_constituents_array()
+
+        if (.not. associated(constituents)) then
+            call endrun('Failed to find variable "constituents"', subname, __LINE__)
+        end if
+
+        call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
+
+        ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index.
+        do i = 1, num_advected
+            ! Vertical index order is reversed between CAM-SIMA and MPAS.
+            do k = 1, pver
+                scalars(i, k, 1:ncells_solve) = &
+                    constituents(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i))
+            end do
+        end do
+
+        nullify(constituents)
+        nullify(scalars)
+
+        ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+        call mpas_dynamical_core % exchange_halo('scalars')
+    end subroutine set_default_constituent
+
+    !> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized
+    !> to prevent physics from attempting to read them from a file. These variables are to be exchanged later
+    !> during dynamics-physics coupling.
+    !> (KCW, 2024-05-23)
+    subroutine mark_variable_as_initialized()
+        character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized'
+        integer :: i
+
+        ! CCPP standard names of physical quantities in the `physics_{state,tend}` derived types.
+        call mark_as_initialized('air_pressure')
+        call mark_as_initialized('air_pressure_at_interface')
+        call mark_as_initialized('air_pressure_of_dry_air')
+        call mark_as_initialized('air_pressure_of_dry_air_at_interface')
+        call mark_as_initialized('air_pressure_thickness')
+        call mark_as_initialized('air_pressure_thickness_of_dry_air')
+        call mark_as_initialized('air_temperature')
+        call mark_as_initialized('dry_static_energy')
+        call mark_as_initialized('eastward_wind')
+        call mark_as_initialized('geopotential_height_wrt_surface')
+        call mark_as_initialized('geopotential_height_wrt_surface_at_interface')
+        call mark_as_initialized('inverse_exner_function_wrt_surface_pressure')
+        call mark_as_initialized('lagrangian_tendency_of_air_pressure')
+        call mark_as_initialized('ln_air_pressure')
+        call mark_as_initialized('ln_air_pressure_at_interface')
+        call mark_as_initialized('ln_air_pressure_of_dry_air')
+        call mark_as_initialized('ln_air_pressure_of_dry_air_at_interface')
+        call mark_as_initialized('northward_wind')
+        call mark_as_initialized('reciprocal_of_air_pressure_thickness')
+        call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air')
+        call mark_as_initialized('surface_air_pressure')
+        call mark_as_initialized('surface_geopotential')
+        call mark_as_initialized('surface_pressure_of_dry_air')
+        call mark_as_initialized('tendency_of_air_temperature_due_to_model_physics')
+        call mark_as_initialized('tendency_of_eastward_wind_due_to_model_physics')
+        call mark_as_initialized('tendency_of_northward_wind_due_to_model_physics')
+
+        ! CCPP standard names of constituents.
+        do i = 1, num_advected
+            call mark_as_initialized(trim(adjustl(const_name(i))))
+        end do
+    end subroutine mark_variable_as_initialized
+
     ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`.
     ! subroutine dyn_run()
     ! end subroutine dyn_run
diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90
index 1effa78c..26b0e345 100644
--- a/src/dynamics/mpas/dyn_grid.F90
+++ b/src/dynamics/mpas/dyn_grid.F90
@@ -1,5 +1,5 @@
 module dyn_grid
-    ! Modules from CAM.
+    ! Modules from CAM-SIMA.
     use cam_abortutils, only: check_allocate, endrun
     use cam_constituents, only: num_advected
     use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, &
@@ -7,8 +7,11 @@ module dyn_grid
                                 max_hcoordname_len
     use cam_initfiles, only: initial_file_get_id
     use cam_map_utils, only: kind_imap => imap
-    use dyn_comp, only: dyn_debug_print, mpas_dynamical_core
-    use dynconst, only: dynconst_init, pi
+    use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, &
+        ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, &
+        ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, &
+        sphere_radius
+    use dynconst, only: constant_pi => pi, rad_to_deg, dynconst_init
     use physics_column_type, only: kind_pcol, physics_column_t
     use physics_grid, only: phys_decomp, phys_grid_init
     use ref_pres, only: ref_pres_init
@@ -26,7 +29,7 @@ module dyn_grid
     implicit none
 
     private
-    ! Provide APIs required by CAM Control.
+    ! Provide APIs required by CAM-SIMA.
     public :: model_grid_init
 
     public :: dyn_grid_id
@@ -40,14 +43,6 @@ module dyn_grid
         'mpas_edge',  &
         'mpas_vertex' &
     ]
-
-    real(kind_r8), parameter :: deg_to_rad = pi / 180.0_kind_r8 ! Convert degrees to radians.
-    real(kind_r8), parameter :: rad_to_deg = 180.0_kind_r8 / pi ! Convert radians to degrees.
-
-    ! Local and global mesh dimensions.
-    integer :: ncells_solve,  nedges_solve,  nvertices_solve
-    integer :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max
-    real(kind_r8) :: sphere_radius
 contains
     !> Initialize various model grids (e.g., dynamics, physics) in terms of dynamics decomposition.
     !> Additionally, MPAS framework initialization and reading time-invariant (e.g., grid/mesh) variables
@@ -57,10 +52,9 @@ module dyn_grid
     ! Called by `cam_init` in `src/control/cam_comp.F90`.
     subroutine model_grid_init()
         character(*), parameter :: subname = 'dyn_grid::model_grid_init'
-        integer, pointer :: ncellssolve => null()
-        integer, pointer :: nedgessolve => null()
-        integer, pointer :: nverticessolve => null()
-        type(file_desc_t), pointer :: pio_file => null()
+        type(file_desc_t), pointer :: pio_file
+
+        nullify(pio_file)
 
         ! Initialize mathematical and physical constants for dynamics.
         call dyn_debug_print('Calling dynconst_init')
@@ -84,20 +78,14 @@ subroutine model_grid_init()
         ! Read time-invariant (e.g., grid/mesh) variables.
         call mpas_dynamical_core % read_write_stream(pio_file, 'r', 'invariant')
 
+        ! Compute local east, north and edge-normal unit vectors whenever time-invariant (e.g., grid/mesh) variables are read.
+        call mpas_dynamical_core % compute_unit_vector()
+
         ! Inquire local and global mesh dimensions and save them as module variables.
         call dyn_debug_print('Inquiring local and global mesh dimensions')
 
-        call mpas_dynamical_core % get_variable_pointer(ncellssolve, 'dim', 'nCellsSolve')
-        call mpas_dynamical_core % get_variable_pointer(nedgessolve, 'dim', 'nEdgesSolve')
-        call mpas_dynamical_core % get_variable_pointer(nverticessolve, 'dim', 'nVerticesSolve')
-
-        ncells_solve = ncellssolve
-        nedges_solve = nedgessolve
-        nvertices_solve = nverticessolve
-
-        nullify(ncellssolve)
-        nullify(nedgessolve)
-        nullify(nverticessolve)
+        call mpas_dynamical_core % get_local_mesh_dimension( &
+            ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels)
 
         call mpas_dynamical_core % get_global_mesh_dimension( &
             ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, &
@@ -138,6 +126,7 @@ end subroutine model_grid_init
     !> Initialize reference pressure by computing necessary variables and calling `ref_pres_init`.
     !> (KCW, 2024-03-25)
     subroutine init_reference_pressure()
+        character(*), parameter :: subname = 'dyn_grid::init_reference_pressure'
         ! Number of pure pressure levels at model top.
         integer, parameter :: num_pure_p_lev = 0
         integer :: ierr
@@ -151,22 +140,24 @@ subroutine init_reference_pressure()
         ! `dzw` denotes the delta/difference between `zw`.
         ! `rdzw` denotes the reciprocal of `dzw`.
         real(kind_r8), allocatable :: zu(:), zw(:), dzw(:)
-        real(kind_r8), pointer :: rdzw(:) => null()
+        real(kind_r8), pointer :: rdzw(:)
+
+        nullify(rdzw)
 
         ! Compute reference height.
         call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw')
 
         allocate(dzw(pver), stat=ierr)
-        call check_allocate(ierr, 'init_reference_pressure', 'dzw', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'dzw(pver)', 'dyn_grid', __LINE__)
 
         dzw(:) = 1.0_kind_r8 / rdzw
 
         nullify(rdzw)
 
         allocate(zw(pverp), stat=ierr)
-        call check_allocate(ierr, 'init_reference_pressure', 'zw', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'zw(pverp)', 'dyn_grid', __LINE__)
         allocate(zu(pver), stat=ierr)
-        call check_allocate(ierr, 'init_reference_pressure', 'zu', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'zu(pver)', 'dyn_grid', __LINE__)
 
         ! In MPAS, zeta coordinates are stored in increasing order (i.e., bottom to top of atmosphere).
         ! In CAM-SIMA, however, index order is reversed (i.e., top to bottom of atmosphere).
@@ -180,12 +171,12 @@ subroutine init_reference_pressure()
 
         ! Compute reference pressure from reference height.
         allocate(p_ref_int(pverp), stat=ierr)
-        call check_allocate(ierr, 'init_reference_pressure', 'p_ref_int', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__)
 
         call std_atm_pres(zw, p_ref_int)
 
         allocate(p_ref_mid(pver), stat=ierr)
-        call check_allocate(ierr, 'init_reference_pressure', 'p_ref_mid', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__)
 
         p_ref_mid(:) = 0.5_kind_r8 * (p_ref_int(1:pver) + p_ref_int(2:pverp))
 
@@ -216,16 +207,22 @@ end subroutine init_reference_pressure
     !> Provide grid and mapping information between global and local indexes to physics by calling `phys_grid_init`.
     !> (KCW, 2024-03-27)
     subroutine init_physics_grid()
+        character(*), parameter :: subname = 'dyn_grid::init_physics_grid'
         character(max_hcoordname_len), allocatable :: dyn_attribute_name(:)
         integer :: hdim1_d, hdim2_d
         integer :: i
         integer :: ierr
-        integer, pointer :: indextocellid(:) => null()  ! Global indexes of cell centers.
-        real(kind_r8), pointer :: areacell(:) => null() ! Cell areas (square meters).
-        real(kind_r8), pointer :: latcell(:) => null()  ! Cell center latitudes (radians).
-        real(kind_r8), pointer :: loncell(:) => null()  ! Cell center longitudes (radians).
+        integer, pointer :: indextocellid(:)  ! Global indexes of cell centers.
+        real(kind_r8), pointer :: areacell(:) ! Cell areas (square meters).
+        real(kind_r8), pointer :: latcell(:)  ! Cell center latitudes (radians).
+        real(kind_r8), pointer :: loncell(:)  ! Cell center longitudes (radians).
         type(physics_column_t), allocatable :: dyn_column(:) ! Grid and mapping information between global and local indexes.
 
+        nullify(areacell)
+        nullify(indextocellid)
+        nullify(latcell)
+        nullify(loncell)
+
         hdim1_d = ncells_global
 
         ! Setting `hdim2_d` to `1` indicates unstructured grid.
@@ -237,7 +234,7 @@ subroutine init_physics_grid()
         call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell')
 
         allocate(dyn_column(ncells_solve), stat=ierr)
-        call check_allocate(ierr, 'init_physics_grid', 'dyn_column', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'dyn_column(ncells_solve)', 'dyn_grid', __LINE__)
 
         do i = 1, ncells_solve
             ! Column information.
@@ -248,7 +245,7 @@ subroutine init_physics_grid()
             ! Cell areas normalized to unit sphere.
             dyn_column(i) % area    = real(areacell(i) / (sphere_radius ** 2), kind_pcol)
             ! Cell weights normalized to unity.
-            dyn_column(i) % weight  = real(areacell(i) / (4.0_kind_r8 * pi * sphere_radius ** 2), kind_pcol)
+            dyn_column(i) % weight  = real(areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2), kind_pcol)
 
             ! File decomposition.
             ! For unstructured grid, `coord_indices` is not used by `phys_grid_init`.
@@ -262,7 +259,8 @@ subroutine init_physics_grid()
             dyn_column(i) % local_dyn_block  = i
             ! `dyn_block_index` is not used due to no dynamics block offset, but it still needs to be allocated.
             allocate(dyn_column(i) % dyn_block_index(0), stat=ierr)
-            call check_allocate(ierr, 'init_physics_grid', 'dyn_column % dyn_block_index', 'dyn_grid', __LINE__)
+            call check_allocate(ierr, subname, 'dyn_column(' // stringify([i]) // ') % dyn_block_index(0)', &
+                'dyn_grid', __LINE__)
         end do
 
         nullify(areacell)
@@ -273,7 +271,7 @@ subroutine init_physics_grid()
         ! `phys_grid_init` expects to receive the `area` attribute from dynamics.
         ! However, do not let it because dynamics grid is different from physics grid.
         allocate(dyn_attribute_name(0), stat=ierr)
-        call check_allocate(ierr, 'init_physics_grid', 'dyn_attribute_name', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'dyn_attribute_name(0)', 'dyn_grid', __LINE__)
 
         call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name)
     end subroutine init_physics_grid
@@ -286,39 +284,50 @@ end subroutine init_physics_grid
     !> * "mpas_vertex": Grid that is centered at MPAS "vertex" points.
     !> (KCW, 2024-03-28)
     subroutine define_cam_grid()
+        character(*), parameter :: subname = 'dyn_grid::define_cam_grid'
         integer :: i
         integer :: ierr
-        integer, pointer :: indextocellid(:) => null()   ! Global indexes of cell centers.
-        integer, pointer :: indextoedgeid(:) => null()   ! Global indexes of edge nodes.
-        integer, pointer :: indextovertexid(:) => null() ! Global indexes of vertex nodes.
-        real(kind_r8), pointer :: areacell(:) => null()  ! Cell areas (square meters).
-        real(kind_r8), pointer :: latcell(:) => null()   ! Cell center latitudes (radians).
-        real(kind_r8), pointer :: latedge(:) => null()   ! Edge node latitudes (radians).
-        real(kind_r8), pointer :: latvertex(:) => null() ! Vertex node latitudes (radians).
-        real(kind_r8), pointer :: loncell(:) => null()   ! Cell center longitudes (radians).
-        real(kind_r8), pointer :: lonedge(:) => null()   ! Edge node longitudes (radians).
-        real(kind_r8), pointer :: lonvertex(:) => null() ! Vertex node longitudes (radians).
+        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_r8), pointer :: areacell(:)  ! Cell areas (square meters).
+        real(kind_r8), pointer :: latcell(:)   ! Cell center latitudes (radians).
+        real(kind_r8), pointer :: latedge(:)   ! Edge node latitudes (radians).
+        real(kind_r8), pointer :: latvertex(:) ! Vertex node latitudes (radians).
+        real(kind_r8), pointer :: loncell(:)   ! Cell center longitudes (radians).
+        real(kind_r8), pointer :: lonedge(:)   ! Edge node longitudes (radians).
+        real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians).
 
         ! Global grid indexes. CAN be safely deallocated because its values are copied into
         ! `cam_grid_attribute_*_t` and `horiz_coord_t`.
         ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
-        integer(kind_imap), pointer :: global_grid_index(:)  => null()
+        integer(kind_imap), pointer :: global_grid_index(:)
         ! Global grid maps. CANNOT be safely deallocated because `cam_filemap_t`
         ! just uses pointers to point at it.
         ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
-        integer(kind_imap), pointer :: global_grid_map(:, :) => null()
+        integer(kind_imap), pointer :: global_grid_map(:, :)
         ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_*_t`
         ! just uses pointers to point at it.
-        real(kind_r8), pointer :: cell_area(:)   => null()
+        real(kind_r8), pointer :: cell_area(:)
         ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_*_t`
         ! just uses pointers to point at it.
-        real(kind_r8), pointer :: cell_weight(:) => null()
+        real(kind_r8), pointer :: cell_weight(:)
         ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
         ! just uses pointers to point at it.
-        type(horiz_coord_t), pointer :: lat_coord => null()
+        type(horiz_coord_t), pointer :: lat_coord
         ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
         ! just uses pointers to point at it.
-        type(horiz_coord_t), pointer :: lon_coord => null()
+        type(horiz_coord_t), pointer :: lon_coord
+
+        nullify(indextocellid, indextoedgeid, indextovertexid)
+        nullify(areacell)
+        nullify(latcell, loncell)
+        nullify(latedge, lonedge)
+        nullify(latvertex, lonvertex)
+
+        nullify(global_grid_index, global_grid_map)
+        nullify(cell_area, cell_weight)
+        nullify(lat_coord, lon_coord)
 
         ! Construct coordinate and grid objects for cell center grid (i.e., "mpas_cell").
         ! Standard MPAS coordinate and dimension names are used.
@@ -329,7 +338,7 @@ subroutine define_cam_grid()
         call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell')
 
         allocate(global_grid_index(ncells_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_index(ncells_solve)', 'dyn_grid', __LINE__)
 
         global_grid_index(:) = int(indextocellid(1:ncells_solve), kind_imap)
 
@@ -339,15 +348,15 @@ subroutine define_cam_grid()
             1, ncells_solve, loncell * rad_to_deg, map=global_grid_index)
 
         allocate(cell_area(ncells_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'cell_area', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'cell_area(ncells_solve)', 'dyn_grid', __LINE__)
         allocate(cell_weight(ncells_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'cell_weight', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'cell_weight(ncells_solve)', 'dyn_grid', __LINE__)
         allocate(global_grid_map(3, ncells_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__)
 
         do i = 1, ncells_solve
             cell_area(i)   = areacell(i)
-            cell_weight(i) = areacell(i) / (4.0_kind_r8 * pi * sphere_radius ** 2)
+            cell_weight(i) = areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2)
 
             global_grid_map(1, i) = int(i, kind_imap)
             global_grid_map(2, i) = int(1, kind_imap)
@@ -363,10 +372,8 @@ subroutine define_cam_grid()
         call cam_grid_attribute_register('mpas_cell', 'cell_weight', 'MPAS cell weight', 'nCells', cell_weight, &
             map=global_grid_index)
 
-        nullify(cell_area)
-        nullify(cell_weight)
-        nullify(lat_coord)
-        nullify(lon_coord)
+        nullify(cell_area, cell_weight)
+        nullify(lat_coord, lon_coord)
 
         ! Construct coordinate and grid objects for cell center grid (i.e., "cam_cell").
         ! Standard CAM-SIMA coordinate and dimension names are used.
@@ -383,14 +390,11 @@ subroutine define_cam_grid()
 
         nullify(areacell)
         nullify(indextocellid)
-        nullify(latcell)
-        nullify(loncell)
+        nullify(latcell, loncell)
 
         deallocate(global_grid_index)
-        nullify(global_grid_index)
-        nullify(global_grid_map)
-        nullify(lat_coord)
-        nullify(lon_coord)
+        nullify(global_grid_index, global_grid_map)
+        nullify(lat_coord, lon_coord)
 
         ! Construct coordinate and grid objects for edge node grid (i.e., "mpas_edge").
         ! Standard MPAS coordinate and dimension names are used.
@@ -400,7 +404,7 @@ subroutine define_cam_grid()
         call mpas_dynamical_core % get_variable_pointer(lonedge, 'mesh', 'lonEdge')
 
         allocate(global_grid_index(nedges_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_index(nedges_solve)', 'dyn_grid', __LINE__)
 
         global_grid_index(:) = int(indextoedgeid(1:nedges_solve), kind_imap)
 
@@ -410,7 +414,7 @@ subroutine define_cam_grid()
             1, nedges_solve, lonedge * rad_to_deg, map=global_grid_index)
 
         allocate(global_grid_map(3, nedges_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__)
 
         do i = 1, nedges_solve
             global_grid_map(1, i) = int(i, kind_imap)
@@ -424,14 +428,11 @@ subroutine define_cam_grid()
             unstruct=.true., block_indexed=.false.)
 
         nullify(indextoedgeid)
-        nullify(latedge)
-        nullify(lonedge)
+        nullify(latedge, lonedge)
 
         deallocate(global_grid_index)
-        nullify(global_grid_index)
-        nullify(global_grid_map)
-        nullify(lat_coord)
-        nullify(lon_coord)
+        nullify(global_grid_index, global_grid_map)
+        nullify(lat_coord, lon_coord)
 
         ! Construct coordinate and grid objects for vertex node grid (i.e., "mpas_vertex").
         ! Standard MPAS coordinate and dimension names are used.
@@ -441,7 +442,7 @@ subroutine define_cam_grid()
         call mpas_dynamical_core % get_variable_pointer(lonvertex, 'mesh', 'lonVertex')
 
         allocate(global_grid_index(nvertices_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_index(nvertices_solve)', 'dyn_grid', __LINE__)
 
         global_grid_index(:) = int(indextovertexid(1:nvertices_solve), kind_imap)
 
@@ -451,7 +452,7 @@ subroutine define_cam_grid()
             1, nvertices_solve, lonvertex * rad_to_deg, map=global_grid_index)
 
         allocate(global_grid_map(3, nvertices_solve), stat=ierr)
-        call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__)
+        call check_allocate(ierr, subname, 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__)
 
         do i = 1, nvertices_solve
            global_grid_map(1, i) = int(i, kind_imap)
@@ -465,14 +466,11 @@ subroutine define_cam_grid()
             unstruct=.true., block_indexed=.false.)
 
         nullify(indextovertexid)
-        nullify(latvertex)
-        nullify(lonvertex)
+        nullify(latvertex, lonvertex)
 
         deallocate(global_grid_index)
-        nullify(global_grid_index)
-        nullify(global_grid_map)
-        nullify(lat_coord)
-        nullify(lon_coord)
+        nullify(global_grid_index, global_grid_map)
+        nullify(lat_coord, lon_coord)
     end subroutine define_cam_grid
 
     !> Helper function for returning grid id given its name.
diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90
index 37a4b120..9d53ca43 100644
--- a/src/dynamics/mpas/stepon.F90
+++ b/src/dynamics/mpas/stepon.F90
@@ -8,7 +8,7 @@ module stepon
     implicit none
 
     private
-    ! Provide APIs required by CAM Control.
+    ! Provide APIs required by CAM-SIMA.
     public :: stepon_init
     public :: stepon_timestep_init
     public :: stepon_run2
diff --git a/src/dynamics/utils/dynconst.F90 b/src/dynamics/utils/dynconst.F90
index c774cf82..78c5b04c 100644
--- a/src/dynamics/utils/dynconst.F90
+++ b/src/dynamics/utils/dynconst.F90
@@ -18,6 +18,11 @@ module dynconst
    !circle's circumference/diameter [unitless]
    real(kind_dyn), parameter, public :: pi = real(phys_pi, kind_dyn)
 
+   ! Convert degrees to radians
+   real(kind_dyn), parameter, public :: deg_to_rad = pi / 180.0_kind_dyn
+   ! Convert radians to degrees
+   real(kind_dyn), parameter, public :: rad_to_deg = 180.0_kind_dyn / pi
+
    ! radius of earth [m]
    real(kind_dyn), protected, public :: rearth
    ! reciprocal of earth's radius [1/m]
@@ -30,6 +35,10 @@ module dynconst
    real(kind_dyn), protected, public :: cpair
    ! Dry air gas constant [J/K/kg]
    real(kind_dyn), protected, public :: rair
+   ! water vapor gas constant [J/K/kg]
+   real(kind_dyn), protected, public :: rh2o
+   ! reference surface pressure [Pa]
+   real(kind_dyn), protected, public :: pref
    ! reference temperature [K]
    real(kind_dyn), protected, public :: tref
    ! reference lapse rate [K/m]
@@ -57,10 +66,12 @@ subroutine dynconst_init
       use physconst, only: phys_omega=>omega
       use physconst, only: phys_cpair=>cpair
       use physconst, only: phys_gravit=>gravit
+      use physconst, only: phys_pref=>pref
       use physconst, only: phys_tref=>tref
       use physconst, only: phys_lapse_rate=>lapse_rate
       use physconst, only: phys_cappa=>cappa
       use physconst, only: phys_rair=>rair
+      use physconst, only: phys_rh2o=>rh2o
 
       !Set constants used by the dynamics:
 
@@ -69,7 +80,9 @@ subroutine dynconst_init
       omega      = real(phys_omega, kind_dyn)
       cpair      = real(phys_cpair, kind_dyn)
       rair       = real(phys_rair, kind_dyn)
+      rh2o       = real(phys_rh2o, kind_dyn)
       gravit     = real(phys_gravit, kind_dyn)
+      pref       = real(phys_pref, kind_dyn)
       tref       = real(phys_tref, kind_dyn)
       lapse_rate = real(phys_lapse_rate, kind_dyn)
       cappa      = real(phys_cappa, kind_dyn)
diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp
index cc42526c..61bd9d3d 160000
--- a/src/physics/ncar_ccpp
+++ b/src/physics/ncar_ccpp
@@ -1 +1 @@
-Subproject commit cc42526c18763b467e3ab67a5c98d8ff2b40bd39
+Subproject commit 61bd9d3dac2abb11bb1e44a2ca34b401da0f44b1