diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam index fe00296b34..ed3c9337f7 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam +++ b/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' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam index 3ff5a2b53f..66c7574789 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam +++ b/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. diff --git a/doc/ChangeLog b/doc/ChangeLog index 07e3f672be..abe108410f 100644 --- a/doc/ChangeLog +++ b/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 diff --git a/src/ionosphere/waccmx/ionosphere_interface.F90 b/src/ionosphere/waccmx/ionosphere_interface.F90 index fa5752f024..59e80260f0 100644 --- a/src/ionosphere/waccmx/ionosphere_interface.F90 +++ b/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 ! diff --git a/src/physics/waccm/mo_aurora.F90 b/src/physics/waccm/mo_aurora.F90 index f6e5039570..7c3443fe4c 100644 --- a/src/physics/waccm/mo_aurora.F90 +++ b/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). !------------------------------------------------------------------------ diff --git a/src/physics/waccmx/steady_state_tei.F90 b/src/physics/waccmx/steady_state_tei.F90 index 135b9055c8..620e314433 100644 --- a/src/physics/waccmx/steady_state_tei.F90 +++ b/src/physics/waccmx/steady_state_tei.F90 @@ -10,30 +10,30 @@ module steady_state_tei use cam_abortutils, only: endrun implicit none - - private + + private !------------------------ - ! PUBLIC: interfaces + ! PUBLIC: interfaces !------------------------ public :: steady_state_tei_init public :: steady_state_tei_tend - + integer :: indxTe = -1 integer :: indxTi = -1 integer :: indxQt = -1 integer :: indxAR = -1 integer :: indxO1 = -1 integer :: indxO2 = -1 - + real(r8) :: op_mass - + contains !============================================================================== subroutine steady_state_tei_init(pbuf2d) - + !----------------------------------------------------------------------- ! Time independent initialization for ionosphere simulation. !----------------------------------------------------------------------- @@ -41,10 +41,8 @@ subroutine steady_state_tei_init(pbuf2d) use constituents, only: cnst_get_ind use mo_chem_utls, only: get_spc_ndx ! Routine to get index of adv_mass array for short lived species use chem_mods, only: adv_mass ! Array holding mass values for short lived species - use infnan, only: nan, assignment(=) type(physics_buffer_desc), pointer :: pbuf2d(:,:) - real(r8) :: nanval call cnst_get_ind( 'O2', indxO2, abort=.true. ) call cnst_get_ind( 'O', indxO1, abort=.true. ) @@ -56,11 +54,7 @@ subroutine steady_state_tei_init(pbuf2d) indxTi = pbuf_get_index( 'TIon' ) indxQt = pbuf_get_index( 'QTeAur' ) indxAR = pbuf_get_index( 'AurIPRateSum' ) - - nanval=nan - call pbuf_set_field(pbuf2d, indxQt, nanval) - call pbuf_set_field(pbuf2d, indxAR, nanval) - + op_mass = adv_mass(get_spc_ndx('Op')) end subroutine steady_state_tei_init @@ -76,7 +70,7 @@ subroutine steady_state_tei_tend(state,istate, dse_tend, pbuf) use perf_mod, only: t_startf, t_stopf ! timing utils use ionos_state_mod, only: ionos_state !------------------------------------------------------------------------------------- - ! Calculate dry static energy and O+ tendency for extended ionosphere simulation + ! Calculate dry static energy and O+ tendency for extended ionosphere simulation !------------------------------------------------------------------------------------- !------------------------------Arguments-------------------------------- @@ -89,29 +83,29 @@ subroutine steady_state_tei_tend(state,istate, dse_tend, pbuf) !---------------------------Local variables------------------------------- real(r8),dimension(pver,state%ncol) :: & - tn, & ! neutral temperature (deg K) - o2, & ! molecular oxygen (mmr) + tn, & ! neutral temperature (deg K) + o2, & ! molecular oxygen (mmr) o1, & ! atomic oxygen (mmr) n2, & ! molecular nitrogen (mmr) ne, & ! electron density (cm3) te, & ! electron temperature (from previous time step) (K) ti, & ! ion temperature (from previous time step) (K) - op, & ! O+ number dens (/cm^3) - o2p, & ! O2+ number dens (/cm^3) - nop, & ! NO+ number dens (/cm^3) + op, & ! O+ number dens (/cm^3) + o2p, & ! O2+ number dens (/cm^3) + nop, & ! NO+ number dens (/cm^3) barm, & ! mean molecular weight (g/mole) qji_ti ! joule heating from qjoule_ti (used ui,vi) ev/g/s real(r8) :: chi(state%ncol) ! solar zenith angle (radians) real(r8) :: qtot(pver,state%ncol) ! total ionization rate (s-1 cm-3) - real(r8) :: pmid(pver,state%ncol) ! mid-level press (dyne/cm^2) ! 10._r8*pmid(:ncol,k) + real(r8) :: pmid(pver,state%ncol) ! mid-level press (dyne/cm^2) ! 10._r8*pmid(:ncol,k) real(r8) :: pint(pverp,state%ncol) ! interface press (dyne/cm^2) - real(r8) :: te_out(pver,state%ncol) ! output electron temperature (deg K) + real(r8) :: te_out(pver,state%ncol) ! output electron temperature (deg K) real(r8) :: ti_out(pver,state%ncol) ! output ion temperature (deg K) real(r8) :: qtotal_out(pver,state%ncol) ! heating rate of neutrals - - integer :: lchnk ! Chunk number + + integer :: lchnk ! Chunk number integer :: ncol ! Number of columns in chunk integer :: i, k @@ -135,7 +129,7 @@ subroutine steady_state_tei_tend(state,istate, dse_tend, pbuf) call pbuf_get_field(pbuf, indxTe, te_ptr) call pbuf_get_field(pbuf, indxTi, ti_ptr) - call pbuf_get_field(pbuf, indxQt, qteaur) + call pbuf_get_field(pbuf, indxQt, qteaur) call pbuf_get_field(pbuf, indxAR, aurIPRateSum) do i =1,ncol @@ -163,7 +157,7 @@ subroutine steady_state_tei_tend(state,istate, dse_tend, pbuf) f107, chi, qtot, qteaur(:ncol), alatm(:ncol,lchnk), & istate%dipmag(:ncol,1), pmid, pint, 1,pver,ncol, & te_out, ti_out, qtotal_out ) - + do i =1,ncol te_ptr(i,pver:1:-1) = te_out(1:pver,i) ti_ptr(i,pver:1:-1) = ti_out(1:pver,i)