Merge pull request #1313 from fvitt/steady_state_TiTe

Fix issue with WACCMX ion / electron temperature solver

Author: fvitt

cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam

@@ -2,7 +2,7 @@ mfilt=1,1,1,1,1,1,1,1,1
 ndens=1,1,1,1,1,1,1,1,1
 nhtfrq=-24,-24,-24,-24,-24,-24,-24,-24,-24
 avgflag_pertape = 'A',  'I',  'I',  'A',  'A', 'I', 'A'
-steady_state_ion_elec_temp=.false.
+steady_state_ion_elec_temp = .true.
 ionos_epotential_amie=.true.
 amienh_files = '$DIN_LOC_ROOT/atm/waccm/amie_data/oct27_31_2003_nh.nc'
 amiesh_files = '$DIN_LOC_ROOT/atm/waccm/amie_data/oct27_31_2003_sh.nc'

cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam

@@ -9,3 +9,4 @@ amiesh_files = '$DIN_LOC_ROOT/atm/waccm/amie_data/oct27_31_2003_sh.nc'
 fincl6 = 'prescr_phihm','prescr_efxm','prescr_kevm','prescr_efxp','prescr_kevp','amie_efx_phys','amie_kev_phys'
 oplus_grid = 144,96
 ionos_npes = 120
+steady_state_ion_elec_temp = .true.

doc/ChangeLog

@@ -1,5 +1,74 @@
 ===============================================================
 
+Tag name: cam6_4_094
+Originator(s): fvitt
+Date: 27 May 2025
+One-line Summary: Fix issue with WACCMX ion and electron temperature solver
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1313
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+ Fix Issue with steady state ion / electron temperature solver #1311
+
+ The auroral module was not overwriting ionization rates initialized to NaN in
+ chunks of physics columns that did not have some columns near the auroral oval
+ which occurred in some configurations (e.g., when SE dycore is used).  Consequently,
+ NaNs were passed to the steady state ion / electron temperature solver.
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: cacraigucar
+
+List all files eliminated: N/A
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+M       cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam
+M       cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam
+ - turn on steady state solver in these tests
+
+M       src/ionosphere/waccmx/ionosphere_interface.F90
+ - add checks for NaN values in ion and electron temperature
+ - improve error messages when NaNs are found
+
+M       src/physics/waccm/mo_aurora.F90
+ - initialize pbuf fields to zero since some chunks of columns not near the
+   auroral oval will not get set to real values
+
+M       src/physics/waccmx/steady_state_tei.F90
+ - remove initialization to NaNs of auroral fields (see above)
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam:
+
+  SMS_D_Ln9.ne16_ne16_mg17.QPX2000.derecho_intel.cam-outfrq9s_amie (Overall: DIFF) details:
+ - expected baseline failure due to setting "steady_state_ion_elec_temp = .true."
+
+  SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL) details:
+ - pre-existing failure due to CLM build-namelist error
+
+derecho/nvhpc/aux_cam: All PASS
+
+izumi/nag/aux_cam: All PASS
+
+izumi/gnu/aux_cam: All PASS
+
+Summarize any changes to answers: bit-for-bit unchanged
+
+===============================================================
+===============================================================
+
 Tag name: cam6_4_093
 Originator(s): eaton
 Date: 19 May 2025
@@ -49,7 +118,7 @@ List all files added and what they do:
 . cime_config/testdefs/testmods_dirs/cam/outfrq9s_bwic/shell_commands
 . cime_config/testdefs/testmods_dirs/cam/outfrq9s_bwic/user_nl_cam
   - add analytic_ic_type='moist_baroclinic_wave_dcmip2016' to the outfrq9s
-    testmods 
+    testmods
 
 . cime_config/testdefs/testmods_dirs/cam/pc7_ne3pg3/user_nl_cam
   namelist settings to test PORT w/ cam7-LT
@@ -102,7 +171,7 @@ cime_config/testdefs/testlist_cam.xml
   - the ERR test was added to make sure the CIME resubmit functionality
     works.  This test moved from izumi to derecho where resubmit is most
     likely to be used.
-. remove  
+. remove
   PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal0
   PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal1
   PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_gnu.cam-ttrac_loadbal3
@@ -112,7 +181,7 @@ cime_config/testdefs/testlist_cam.xml
   - physics load balancing only works with the FV dycore.  These tests
     aren't doing anything.
 . add comments to simple model tests
-  
+
 
 src/physics/cam/physpkg.F90
 src/physics/cam7/physpkg.F90

src/ionosphere/waccmx/ionosphere_interface.F90

@@ -840,6 +840,11 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
             end do
          end do
 
+         if (state_debug_checks) then
+            call shr_assert_in_domain(te_blck, is_nan=.false., varname="te_blck", msg="NaN found in te_blck in ionosphere_run2")
+            call shr_assert_in_domain(ti_blck, is_nan=.false., varname="ti_blck", msg="NaN found in ti_blck in ionosphere_run2")
+         end if
+
          call t_startf('d_pie_coupling')
 
          ! Compute geometric height and some diagnostic fields needed by
@@ -861,10 +866,10 @@ subroutine ionosphere_run2(phys_state, pbuf2d)
          call t_stopf ('d_pie_coupling')
 
          if (state_debug_checks) then
-            call shr_assert_in_domain(ui_blck, is_nan=.false., varname="ui_blck", msg="NaN found in ionosphere_run2")
-            call shr_assert_in_domain(vi_blck, is_nan=.false., varname="vi_blck", msg="NaN found in ionosphere_run2")
-            call shr_assert_in_domain(wi_blck, is_nan=.false., varname="wi_blck", msg="NaN found in ionosphere_run2")
-            call shr_assert_in_domain(opmmr_blck, is_nan=.false., varname="opmmr_blck", msg="NaN found in ionosphere_run2")
+            call shr_assert_in_domain(ui_blck, is_nan=.false., varname="ui_blck", msg="NaN found in ui_blck in ionosphere_run2")
+            call shr_assert_in_domain(vi_blck, is_nan=.false., varname="vi_blck", msg="NaN found in vi_blck in ionosphere_run2")
+            call shr_assert_in_domain(wi_blck, is_nan=.false., varname="wi_blck", msg="NaN found in wi_blck in ionosphere_run2")
+            call shr_assert_in_domain(opmmr_blck, is_nan=.false., varname="opmmr_blck", msg="NaN found in opmmr_blck in ionosphere_run2")
          end if
 
          !

src/physics/waccm/mo_aurora.F90

@@ -5,7 +5,7 @@ module mo_aurora
 ! Auroral oval parameterization. See reference:
 ! R.G. Roble, E.C. Ridley
 ! An auroral model for the NCAR thermospheric general circulation model (TGCM)
-! Annales Geophysicae,5A, (6), 369-382, 1987. 
+! Annales Geophysicae,5A, (6), 369-382, 1987.
 !
 ! The aurora oval is a circle in auroral circle coordinates.  Auroral circle
 !  coordinates are offset from magnetic coordinates by offa degrees (radians)
@@ -46,7 +46,7 @@ module mo_aurora
 !   1) sub aurora_cons called once per time step from advance.
 !   2) sub aurora called from dynamics, inside parallel latitude scan.
 !   3) subs aurora_cusp and aurora_heat called from sub aurora.
-!   4) sub aurora_ions called from sub aurora. 
+!   4) sub aurora_ions called from sub aurora.
 !
 !-----------------------------------------------------------------------
 
@@ -71,7 +71,7 @@ module mo_aurora
       private
       public :: aurora_inti, aurora_timestep_init, aurora
       public :: aurora_register
-      
+
       integer, parameter  :: isouth = 1
       integer, parameter  :: inorth = 2
 
@@ -143,7 +143,6 @@ module mo_aurora
 
       contains
 
-        
       !----------------------------------------------------------------------
       !----------------------------------------------------------------------
       subroutine aurora_register
@@ -197,9 +196,14 @@ subroutine aurora_inti(pbuf2d)
          call pbuf_set_field(pbuf2d, indxKev, x_nan)
       endif
 
+      ! initialize these fields to zero here since some chunks of colunms not near
+      ! the poles will not get set to real values
       if (indxAIPRS>0) then
          call pbuf_set_field(pbuf2d, indxAIPRS, 0._r8)
       endif
+      if (indxQTe>0) then
+         call pbuf_set_field(pbuf2d, indxQTe, 0._r8)
+      endif
 
       theta0(:) = nan
       offa(:) = nan
@@ -347,9 +351,9 @@ subroutine aurora_timestep_init( )
 !-----------------------------------------------------------------------
 
       rh = (h2 - h1) / (h1 + h2)
-      h0     = 0.5_r8 * (h1 + h2) * d2r   
-      
-      
+      h0     = 0.5_r8 * (h1 + h2) * d2r
+
+
 ! roth = MLT of max width of aurora in hours
 ! rote = MLT of max energy flux of aurora in hours
 
@@ -459,7 +463,6 @@ subroutine aurora_prod( tn, o2, o1, mbar, rlats, &
       logical  :: do_aurora(ncol)
 
       real(r8) :: dayfrac, rotation
-      real(r8), pointer :: qteaur(:)    ! for electron temperature
 
       if (.not. aurora_active) return
 
@@ -479,11 +482,6 @@ subroutine aurora_prod( tn, o2, o1, mbar, rlats, &
       call outfld( 'ALONM', r2d*alonm(:ncol,lchnk), ncol, lchnk )
       call outfld( 'ALATM', r2d*alatm(:ncol,lchnk), ncol, lchnk )
 
-      if (indxQTe>0) then
-        call pbuf_get_field(pbuf, indxQTe, qteaur)
-        qteaur(:) = 0._r8
-     endif
-     
 !-----------------------------------------------------------------------
 !    aurora is active for columns poleward of 30 deg
 !-----------------------------------------------------------------------
@@ -501,7 +499,7 @@ subroutine aurora_prod( tn, o2, o1, mbar, rlats, &
       do i = 1,ncol
         if( do_aurora(i) ) then
           dlat_aur(i) = alatm(i,lchnk)
-          dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it 
+          dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it
           if( dlon_aur(i) > pi ) then
              dlon_aur(i) = dlon_aur(i) - twopi
           else if( dlon_aur(i) < -pi ) then
@@ -659,7 +657,7 @@ subroutine aurora_hrate( tn, mbar, rlats, &
       do i = 1,ncol
         if( do_aurora(i) ) then
           dlat_aur(i) = alatm(i,lchnk)
-          dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it 
+          dlon_aur(i) = alonm(i,lchnk) + rotation ! rotate it
           if( dlon_aur(i) > pi ) then
              dlon_aur(i) = dlon_aur(i) - twopi
           else if( dlon_aur(i) < -pi ) then
@@ -714,7 +712,7 @@ subroutine aurora_hrate( tn, mbar, rlats, &
       call aurora_heat( flux, flux2, alfa, alfa2, &
                         drizl, do_aurora, hemis, &
                         alon, colat, ncol, pbuf )
- 
+
 !-----------------------------------------------------------------------
 ! 	... auroral additions to ionization rates
 !-----------------------------------------------------------------------
@@ -763,7 +761,7 @@ subroutine aurora_cusp( cusp, do_aurora, hemis, colat, alon, ncol )
          cusp(:) = 0._r8
       endwhere
 
-      end subroutine aurora_cusp 
+      end subroutine aurora_cusp
 
       subroutine aurora_heat( flux, flux2, alfa, alfa2, &
                               drizl, do_aurora, hemis, &
@@ -772,7 +770,7 @@ subroutine aurora_heat( flux, flux2, alfa, alfa2, &
 ! 	... calculate alfa, flux, and drizzle
 !-----------------------------------------------------------------------
       use physics_buffer,only: physics_buffer_desc,pbuf_get_field
-      
+
       implicit none
 
 !-----------------------------------------------------------------------
@@ -798,12 +796,12 @@ subroutine aurora_heat( flux, flux2, alfa, alfa2, &
         halfwidth, &                            ! oval half-width
         wrk, &                                  ! temp wrk array
         dtheta                                  ! latitudinal variation (Gaussian)
-      real(r8) :: ekev 
+      real(r8) :: ekev
       real(r8), pointer   :: pr_efx(:) ! Pointer to pbuf prescribed energy flux (mW m-2)
       real(r8), pointer   :: pr_kev(:) ! Pointer to pbuf prescribed mean energy (keV)
       real(r8), pointer   :: qteaur(:)    ! for electron temperature
       integer  :: n
-      
+
 !-----------------------------------------------------------------------
 ! Low-energy protons:
 !
@@ -850,7 +848,7 @@ subroutine aurora_heat( flux, flux2, alfa, alfa2, &
       endwhere
 
 !-----------------------------------------------------------------------
-! 	... for electron temperature (used in settei):  
+! 	... for electron temperature (used in settei):
 !-----------------------------------------------------------------------
       if (indxQTe>0) then
          call pbuf_get_field(pbuf, indxQTe, qteaur)
@@ -954,14 +952,14 @@ subroutine aurora_ions( drizl, cusp, alfa1, alfa2, &
       real(r8) :: wrk(ncol,pver)
 
       real(r8), pointer   :: aurIPRateSum(:,:) ! Pointer to pbuf auroral ion production sum for O2+,O+,N2+ (s-1 cm-3)
-      
+
       qia(:) = 0._r8
       wrk(:,:) = 0._r8
 
       !-----------------------------------------------------------
-      !  Point to production rates array in physics buffer where 
-      !  rates will be stored for ionosphere module access.  Also, 
-      !  initialize rates to zero before column loop since only 
+      !  Point to production rates array in physics buffer where
+      !  rates will be stored for ionosphere module access.  Also,
+      !  initialize rates to zero before column loop since only
       !  daylight values are filled
       !-----------------------------------------------------------
       if (indxAIPRS>0) then
@@ -1041,15 +1039,15 @@ subroutine aurora_ions( drizl, cusp, alfa1, alfa2, &
       end do level_loop
 
       !----------------------------------------------------------------
-      !  Store the sum of the ion production rates in pbuf to be used 
-      !  in the ionosx module 
+      !  Store the sum of the ion production rates in pbuf to be used
+      !  in the ionosx module
       !----------------------------------------------------------------
       if (indxAIPRS>0) then
-      
-        aurIPRateSum(1:ncol,1:pver) = wrk(1:ncol,1:pver) 
-      
+
+        aurIPRateSum(1:ncol,1:pver) = wrk(1:ncol,1:pver)
+
       endif
-  
+
       call outfld( 'QSUM', wrk, ncol, lchnk )
 
       end subroutine aurora_ions
@@ -1165,9 +1163,9 @@ subroutine aion( si, so, do_aurora, ncol )
 !-----------------------------------------------------------------------
 ! Calculates integrated f(x) needed for total auroral ionization.
 ! See equations (10-12) in Roble,1987.
-! Coefficients for equation (12) of Roble,1987 are in variable cc 
+! Coefficients for equation (12) of Roble,1987 are in variable cc
 ! (revised since 1987):
-! Uses the identity x**y = exp(y*ln(x)) for performance 
+! Uses the identity x**y = exp(y*ln(x)) for performance
 ! (fewer (1/2) trancendental functions are required).
 !------------------------------------------------------------------------