Skip to content

cam6_4_094: Fix issue with WACCMX ion / electron temperature solver #1313

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
May 27, 2025
Merged

cam6_4_094: Fix issue with WACCMX ion / electron temperature solver #1313

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
30a75f6
fix issue with NaNs in AurIPRateSum
fvitt May 21, 2025
d45ec20
ChangeLog draft
fvitt May 27, 2025
dd5f7e6
ChangeLog update
fvitt May 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion cime_config/testdefs/testmods_dirs/cam/outfrq1d_amie/user_nl_cam
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down
1 change: 1 addition & 0 deletions cime_config/testdefs/testmods_dirs/cam/outfrq9s_amie/user_nl_cam
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
13 changes: 9 additions & 4 deletions src/ionosphere/waccmx/ionosphere_interface.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

!
Expand Down
64 changes: 31 additions & 33 deletions src/physics/waccm/mo_aurora.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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.
!
!-----------------------------------------------------------------------

Expand All @@ -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

Expand Down Expand Up @@ -143,7 +143,6 @@ module mo_aurora

contains


!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine aurora_register
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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
!-----------------------------------------------------------------------
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
!-----------------------------------------------------------------------
Expand Down Expand Up @@ -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, &
Expand All @@ -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

!-----------------------------------------------------------------------
Expand All @@ -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:
!
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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).
!------------------------------------------------------------------------

Expand Down
46 changes: 20 additions & 26 deletions src/physics/waccmx/steady_state_tei.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,41 +10,39 @@ 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.
!-----------------------------------------------------------------------

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. )
Expand All @@ -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
Expand All @@ -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--------------------------------
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down