From 4f35eebb66e17a62d2793ba600d8b2f2a968ee78 Mon Sep 17 00:00:00 2001
From: Bill Skamarock <skamaroc@ucar.edu>
Date: Fri, 17 Oct 2025 13:30:56 -0600
Subject: [PATCH] This commit provides the capability of turning off use of the
 microphysics tendency above a user-specified height.  This change alleviates
 some instabilities encountered in deep-domain MPAS simulations (model tops
 above the stratopause).

The changes include:

(1) addition of a namelist variable specifying the height above which the
microphysics tendencies are to be ignored, and an integer variable that is set
when the integration starts specifying the level based on that height.
(2) additions to src/core_atmosphere/physics/mpas_atmphys_init.F to set the
level based on the namelist-specified height.
(3) changes in src/core_atmosphere/physics/mpas_atmphys_interface.F that disable
the microphysics tendancy updates above the specified level.
---
 src/core_atmosphere/Registry.xml              |  8 ++
 .../physics/mpas_atmphys_init.F               | 20 +++-
 .../physics/mpas_atmphys_interface.F          | 96 ++++++++++++++-----
 3 files changed, 98 insertions(+), 26 deletions(-)

diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 4281c40bba..b0bd5e9e31 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -1599,6 +1599,9 @@
                 <var name="nearestRelaxationCell" type="integer" dimensions="nCells" units="-"
                      description="For cells in the specified zone, gives the index of the nearest cell in the relaxation zone"/>
 
+                <var name="mp_top_level" type="integer" dimensions="" units="-"
+                     description="Vertical layer above which microphysics tendency is ignored"/>
+
         </var_struct>
 
         <var_struct name="state" time_levs="2">
@@ -2215,6 +2218,11 @@
                      description="number of microphysics time-steps per physics time-steps"
                      possible_values="Positive integers"/>
 
+                <nml_option name="config_microphysics_top" type="real" default_value="45000." in_defaults="false"
+                     units="m"
+                     description="height above which to exclude microphysics tendencies"
+                     possible_values="Positive reals"/>
+
                 <nml_option name="config_radtlw_interval" type="character" default_value="00:30:00"
                      units="-"
                      description="time interval between calls to parameterization of long-wave radiation"
diff --git a/src/core_atmosphere/physics/mpas_atmphys_init.F b/src/core_atmosphere/physics/mpas_atmphys_init.F
index 5e851ad923..32969cb5df 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_init.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_init.F
@@ -73,7 +73,8 @@ module mpas_atmphys_init
 !   Laura D. Fowler (laura@ucar.edu) / 2024-02-14.
 ! * added call to subroutine init_lsm_noahmp to initialize the Noah-MP land surface scheme.
 !   Laura D. Fowler (laura@ucar.edu) / 2024-03-11.
-
+! * added initialization of the integer variable mp_top_level to the physics_init code.
+!   Bill Skamarock / 2025-10-17 
 
  contains
 
@@ -149,13 +150,16 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
  integer,dimension(:),pointer:: jindx1_tau, jindx2_tau
  real(kind=RKIND),dimension(:),pointer:: ddy_j1tau, ddy_j2tau
  real(kind=RKIND),dimension(:),pointer:: latCell
+ real(kind=RKIND),pointer:: config_microphysics_top
+ integer,pointer:: mp_top_level
 
 !local variables and arrays:
  type(MPAS_Time_Type):: currTime
 
  logical:: init_done
  integer:: ierr,julday 
- integer:: iCell,iLag
+ integer:: iCell,iLag,k
+ real(kind=RKIND):: layer_height
 
 !-----------------------------------------------------------------------------------------------------------------
 !call mpas_log_write('')
@@ -170,6 +174,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
  call mpas_pool_get_config(configs,'config_sfclayer_scheme'  ,config_sfclayer_scheme  )
  call mpas_pool_get_config(configs,'config_radt_lw_scheme'   ,config_radt_lw_scheme   )
  call mpas_pool_get_config(configs,'config_radt_sw_scheme'   ,config_radt_sw_scheme   )
+ call mpas_pool_get_config(configs,'config_microphysics_top' ,config_microphysics_top )
 
  call mpas_pool_get_dimension(mesh,'nCellsSolve',nCellsSolve)
  call mpas_pool_get_dimension(mesh,'nLags'      ,nLags      )
@@ -237,6 +242,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
  call mpas_pool_get_dimension(mesh,'nVertLevels',nVertLevels  )
  call mpas_pool_get_array(mesh,'rdzw'   ,rdzw   )
  call mpas_pool_get_array(mesh,'dzu'    ,dzu    )
+ call mpas_pool_get_array(mesh,'mp_top_level'    ,mp_top_level    )
 
  currTime = mpas_get_clock_time(clock,MPAS_NOW,ierr)
  call mpas_get_time(curr_time=currTime,DoY=julday,ierr=ierr)
@@ -265,6 +271,16 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
     enddo
  endif
 
+! set the level above which microphysics tendencies will be ignored
+ mp_top_level = 0
+ layer_height = 0.5/rdzw(1)
+ if(layer_height .le. config_microphysics_top) mp_top_level = 1
+ do k=2,nVertLevels
+    layer_height = layer_height + dzu(k)
+    if(layer_height .le. config_microphysics_top) mp_top_level = k
+ end do
+ call mpas_log_write('physics init: config_microphysics_top and level $r $i ',realArgs=(/config_microphysics_top/),intArgs=(/mp_top_level/)) 
+
 !initialization of counters i_acsw* and i_aclw*. i_acsw* and i_aclw* track the number of times
 !the accumulated long and short-wave radiation fluxes exceed their prescribed theshold values.
  if(.not. config_do_restart) then
diff --git a/src/core_atmosphere/physics/mpas_atmphys_interface.F b/src/core_atmosphere/physics/mpas_atmphys_interface.F
index 71e46dfcd2..58206cf9c9 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_interface.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_interface.F
@@ -64,7 +64,8 @@ module mpas_atmphys_interface
 !   Laura D. Fowler (laura@ucar.edu) / 2016-04-11.
 ! * corrected the calculation of the surface pressure, mainly extrapolation of the air density to the surface.
 !   Laura D. Fowler (laura@ucar.edu) / 2016-04-25.
-
+! * changed microphysics update such that the microphsyics tendency is ignored above level mp_top_level
+!   Bill Skamarock / 2025/10/17 
 
  contains
 
@@ -810,6 +811,7 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
  real(kind=RKIND),dimension(:,:),pointer  :: rthmpten,rqvmpten,rqcmpten,rqrmpten,rqimpten,rqsmpten,rqgmpten
  real(kind=RKIND),dimension(:,:),pointer  :: rncmpten,rnimpten,rnrmpten,rnifampten,rnwfampten
  real(kind=RKIND),dimension(:,:,:),pointer:: scalars
+ integer,pointer:: mp_top_level
 
 !local variables:
  integer:: icount
@@ -822,6 +824,7 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
 
  call mpas_pool_get_array(mesh,'zz'   ,zz   )
  call mpas_pool_get_array(mesh,'zgrid',zgrid)
+ call mpas_pool_get_array(mesh,'mp_top_level',mp_top_level)
 
  call mpas_pool_get_array(diag,'exner'           ,exner           )
  call mpas_pool_get_array(diag,'exner_base'      ,exner_b         )
@@ -849,7 +852,7 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
 
 !update variables needed in the dynamical core:
  do j = jts,jte
- do k = kts,kte
+ do k = kts,mp_top_level
  do i = its,ite
 
     !initializes tendency of coupled potential temperature potential temperature, and
@@ -859,11 +862,11 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
 
     !updates water vapor, cloud liquid water, rain mixing ratios, modified potential temperature,
     !tendency of coupled potential temperature, and potential temperature heating rate from microphysics:
+
     qv(k,i) = qv_p(i,k,j)
+    theta_m(k,i) = th_p(i,k,j) * (1._RKIND+rvord*qv_p(i,k,j))
     qc(k,i) = qc_p(i,k,j)
     qr(k,i) = qr_p(i,k,j)
-
-    theta_m(k,i) = th_p(i,k,j) * (1._RKIND+rvord*qv_p(i,k,j))
     rt_diabatic_tend(k,i) = (theta_m(k,i) - rt_diabatic_tend(k,i))/dt_dyn
     dtheta_dt_mp(k,i)     = (theta_m(k,i)/(1._RKIND+rvord*qv(k,i))-dtheta_dt_mp(k,i))/(dt_dyn)
 
@@ -881,7 +884,20 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
  enddo
  enddo
 
-!update surface pressure and calculates the surface pressure tendency:
+ ! no microphysics updates above mp_top_level
+ ! why are we not applying this to theta and qv?
+ 
+ if(kte .gt. mp_top_level) then
+    do i = its,ite
+       do k = mp_top_level+1,kte
+          qc(k,i) = 0.
+          qr(k,i) = 0.
+          rt_diabatic_tend(k,i) = 0.
+       end do
+    enddo
+ endif
+ 
+ ! update surface pressure and calculates the surface pressure tendency:
  do j = jts,jte
  do i = its,ite
     tem1 = zgrid(2,i)-zgrid(1,i)
@@ -917,22 +933,38 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
        call mpas_pool_get_array(diag_physics,'refl10cm',refl10cm)
 
        do j = jts,jte
-       do k = kts,kte
+       do k = kts,mp_top_level
        do i = its,ite
-          qi(k,i) = qi_p(i,k,j)
-          qs(k,i) = qs_p(i,k,j)
-          qg(k,i) = qg_p(i,k,j)
-
-          rainprod(k,i) = rainprod_p(i,k,j)
-          evapprod(k,i) = evapprod_p(i,k,j)
-          re_cloud(k,i) = recloud_p(i,k,j)
-          re_ice(k,i)   = reice_p(i,k,j)
-          re_snow(k,i)  = resnow_p(i,k,j)
+           qi(k,i) = qi_p(i,k,j)
+           qs(k,i) = qs_p(i,k,j)
+           qg(k,i) = qg_p(i,k,j)
+           rainprod(k,i) = rainprod_p(i,k,j)
+           evapprod(k,i) = evapprod_p(i,k,j)
+           re_cloud(k,i) = recloud_p(i,k,j)
+           re_ice(k,i)   = reice_p(i,k,j)
+           re_snow(k,i)  = resnow_p(i,k,j)
        enddo
        enddo
        enddo
 
-       mp2_select: select case(trim(mp_scheme))
+       ! no microphysics updates above mp_top_level
+       
+       if(kte .gt. mp_top_level) then
+          do i = its,ite
+             do k = mp_top_level+1,kte
+                qi(k,i) = 0.
+                qs(k,i) = 0.
+                qg(k,i) = 0.
+                rainprod(k,i) = 0.
+                evapprod(k,i) = 0.
+                re_cloud(k,i) = re_qc_bg
+                re_ice(k,i)   = re_qi_bg
+                re_snow(k,i)  = re_qs_bg
+             end do
+          end do
+       end if
+
+         mp2_select: select case(trim(mp_scheme))
           case("mp_thompson","mp_thompson_aerosols")
              call mpas_pool_get_dimension(state,'index_ni',index_ni)
              call mpas_pool_get_dimension(state,'index_nr',index_nr)
@@ -993,18 +1025,34 @@ subroutine microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,te
        call mpas_pool_get_array(tend_physics,'rqsmpten',rqsmpten)
        call mpas_pool_get_array(tend_physics,'rqgmpten',rqgmpten)
 
-       do k = kts,kte
        do i = its,ite
-          rthmpten(k,i) = (theta_m(k,i)/(1._RKIND+R_v/R_d*max(0._RKIND,qv(k,i)))-rthmpten(k,i))/dt_dyn
-          rqvmpten(k,i) = (qv(k,i)-rqvmpten(k,i))/dt_dyn
-          rqcmpten(k,i) = (qc(k,i)-rqcmpten(k,i))/dt_dyn
-          rqrmpten(k,i) = (qr(k,i)-rqrmpten(k,i))/dt_dyn
-          rqimpten(k,i) = (qi(k,i)-rqimpten(k,i))/dt_dyn
-          rqsmpten(k,i) = (qs(k,i)-rqsmpten(k,i))/dt_dyn
-          rqgmpten(k,i) = (qg(k,i)-rqgmpten(k,i))/dt_dyn
+       do k = kts,mp_top_level
+           rthmpten(k,i) = (theta_m(k,i)/(1._RKIND+R_v/R_d*max(0._RKIND,qv(k,i)))-rthmpten(k,i))/dt_dyn
+           rqvmpten(k,i) = (qv(k,i)-rqvmpten(k,i))/dt_dyn
+           rqcmpten(k,i) = (qc(k,i)-rqcmpten(k,i))/dt_dyn
+           rqrmpten(k,i) = (qr(k,i)-rqrmpten(k,i))/dt_dyn
+           rqimpten(k,i) = (qi(k,i)-rqimpten(k,i))/dt_dyn
+           rqsmpten(k,i) = (qs(k,i)-rqsmpten(k,i))/dt_dyn
+           rqgmpten(k,i) = (qg(k,i)-rqgmpten(k,i))/dt_dyn
+!changes from Wei
        enddo
        enddo
 
+       ! no microphysics tendency above level mp_top_level
+       if(kte .gt. mp_top_level) then
+          do i = its, ite
+             do k = mp_top_level+1, kte
+                rthmpten(k,i) = 0.
+                rqvmpten(k,i) = 0.
+                rqcmpten(k,i) = 0.
+                rqrmpten(k,i) = 0.
+                rqimpten(k,i) = 0.
+                rqsmpten(k,i) = 0.
+                rqgmpten(k,i) = 0.
+             end do
+          end do
+       endif
+
        mp2_tend_select: select case(trim(mp_scheme))
           case("mp_thompson","mp_thompson_aerosols")
              call mpas_pool_get_array(tend_physics,'rnimpten',rnimpten)