Skip to content

Remove increased divergence damping for SE dycore + frontogenesis function bug fix for FV dycore #1455

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

Closed
PeterHjortLauritzen wants to merge 8 commits into from
Closed

Remove increased divergence damping for SE dycore + frontogenesis function bug fix for FV dycore #1455

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
48a9b71
fix bug in diagnostics in prim_state
PeterHjortLauritzen Dec 15, 2025
e85ec68
remove increased divergence damping - science change
PeterHjortLauritzen Dec 15, 2025
2d6e1a3
update ChangeLog
PeterHjortLauritzen Dec 15, 2025
8a5f2c9
changes from Jack Chen for bug in frontogenesis code with FV
PeterHjortLauritzen Dec 15, 2025
72ca9e6
remove old 70 level WACCM6 SE dycore test
PeterHjortLauritzen Dec 15, 2025
3bcc268
update ChangeLog
PeterHjortLauritzen Dec 15, 2025
822da0e
more ChangeLog changes
PeterHjortLauritzen Dec 15, 2025
af2e2c3
minor diagnostics mod
PeterHjortLauritzen Dec 16, 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
17 changes: 8 additions & 9 deletions bld/namelist_files/namelist_defaults_cam.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3168,21 +3168,19 @@

<se_hypervis_scaling se_refined_mesh="1" hypervis_type="tensor" >3.22D0 </se_hypervis_scaling>

<se_hypervis_subcycle > 3 </se_hypervis_subcycle>
<se_hypervis_subcycle > 1 </se_hypervis_subcycle>
<se_hypervis_subcycle waccm_phys="1" > 2 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" > 4 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" waccm_phys="1" model_top="none"> 10 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" waccm_phys="1" waccmx="0" model_top="none"> 9 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" waccm_phys="1" waccmx="0" model_top="none"> 8 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" waccm_phys="1" waccmx="1" model_top="none"> 5 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="lt" > 2 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="mt" > 3 </se_hypervis_subcycle>
<se_hypervis_subcycle se_refined_mesh="1" > 3 </se_hypervis_subcycle>
<se_hypervis_subcycle se_refined_mesh="1" > 2 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne0np4CONUS.ne30x8" waccm_phys="1" > 1 </se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" model_top="ht" > 18</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" model_top="xt" > 18</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="ht" > 18</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="xt" > 18</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" model_top="ht" > 6</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne16np4" model_top="xt" >18</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="ht" > 6</se_hypervis_subcycle>
<se_hypervis_subcycle hgrid="ne30np4" model_top="xt" >18</se_hypervis_subcycle>

<se_hypervis_subcycle_sponge > 1 </se_hypervis_subcycle_sponge>
<se_hypervis_subcycle_sponge hgrid="ne30np4" model_top="mt" > 3 </se_hypervis_subcycle_sponge>
Expand All @@ -3198,7 +3196,6 @@
<se_hypervis_subcycle_sponge hgrid="ne30np4" model_top="xt" > 40 </se_hypervis_subcycle_sponge>

<se_hypervis_subcycle_q>1 </se_hypervis_subcycle_q>
<se_hypervis_subcycle_q hgrid="ne16np4">2 </se_hypervis_subcycle_q>

<se_limiter_option> 8 </se_limiter_option>

Expand All @@ -3219,6 +3216,8 @@
<se_nu_div waccmx="1"> 10.e15 </se_nu_div>

<se_nu_p>-1 </se_nu_p>
<se_nu_p hgrid="ne16np4" waccm_phys="1" waccmx="0">6.e15</se_nu_p>
<se_nu_p hgrid="ne16np4" model_top="ht" >6.e15</se_nu_p>

<se_nu_top > 1.25e5 </se_nu_top>
<se_nu_top waccmx="1" > 1.0e6 </se_nu_top>
Expand Down
9 changes: 0 additions & 9 deletions cime_config/testdefs/testlist_cam.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2605,15 +2605,6 @@
<option name="comment" >WACCM7 with MA chemistry</option>
</options>
</test>
<test compset="FW2000climo" grid="ne30pg3_ne30pg3_mg17" name="SMS_Ln9" testmods="cam/outfrq9s_rrtmgp">
<machines>
<machine name="derecho" compiler="intel" category="aux_cam"/>
</machines>
<options>
<option name="wallclock">00:20:00</option>
<option name="comment">WACCM w/ RRTMGP</option>
</options>
</test>
<test compset="FW2000climo" grid="ne30pg3_ne30pg3_mg17" name="SMS_D_Ln9" testmods="cam/outfrq9s">
<machines>
<machine name="derecho" compiler="intel" category="waccm"/>
Expand Down
35 changes: 35 additions & 0 deletions doc/ChangeLog
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,41 @@

===============================================================

Tag name: cam6_4_???
Originator(s): pel,cchen
Date:
One-line Summary: remove increased divergence damping in SE dycore + frontogenesis fct for FV
Github PR URL: https://github.com/ESCOMP/CAM/pull/

Purpose of changes (include the issue number and title text for each relevant GitHub issue):
remove increased divergence damping in SE dycore: https://github.com/ESCOMP/CAM/issues/1449
fix bug in frontogenesis function for FV dycore: https://github.com/ESCOMP/CAM/issues/1414
(science tested by ACOM; there will be namelist changes to re-tune QBO with the fixed frontogenesis function)
remove failing WACCM6 70 level CAM6 test with SE dycore:
SMS_Ln9.ne30pg3_ne30pg3_mg17.FW2000climo.derecho_intel.cam-outfrq9s_rrtmgp
https://github.com/ESCOMP/CAM/issues/1450

All SE regression tests will have answer changes
All FV WACCM regression test will have answer changes

Describe any changes made to build system: N/A

Describe any changes made to the namelist: defaults for se_nu_p, se_hypervis_subcycle, ... changed

List any changes to the defaults for the boundary datasets: N/A

Describe any substantial timing or memory changes: N/A

Code reviewed by:

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:

===============================================================

Tag name: cam6_4_089
Originator(s): cacraig, jedwards
Date: April 28, 2025
Expand Down
95 changes: 95 additions & 0 deletions src/dynamics/fv/gravity_waves_sources.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,11 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
real(r8) :: ug(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! U-wind ghosted
real(r8) :: vg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! V-wind ghosted

real(r8) :: pmg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! mid-point pressure ghosted
real(r8) :: dptdp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)
real(r8) :: dudp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)
real(r8) :: dvdp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)

real(r8) :: tglat ! tangent-latitude
integer :: i,j,k
integer :: im, ip
Expand Down Expand Up @@ -100,6 +105,12 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
call ghost_array(grid, u3, ug)
call ghost_array(grid, v3, vg)

call ghost_array(grid, pm, pmg)

call compute_vertical_derivative(grid,pe,pm,ptemp,dptdp)
call compute_vertical_derivative(grid,pe,pm,u3,dudp)
call compute_vertical_derivative(grid,pe,pm,v3,dvdp)

!$omp parallel do private (i,j,k)
do k=1, plev
do i=beglonxy, endlonxy
Expand All @@ -109,14 +120,26 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
pty(i,k,j) = ( ptg(i,k,j+1) - ptg(i,k,j-1) ) / (2._r8 * grid%dp)
pty(i,k,j) = pty(i,k,j) / aearth

! Topography correction term
pty(i,k,j) = pty(i,k,j) - dptdp(i,k,j) * &
( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth

! U-wind
uy(i,k,j) = ( ug(i,k,j+1) - ug(i,k,j-1) ) / (2._r8 * grid%dp)
uy(i,k,j) = uy(i,k,j) / aearth

! Topography correction term
uy(i,k,j) = uy(i,k,j) - dudp(i,k,j) * &
( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth

! V-wind
vy(i,k,j) = ( vg(i,k,j+1) - vg(i,k,j-1) ) / (2._r8 * grid%dp)
vy(i,k,j) = vy(i,k,j) / aearth

! Topography correction term
vy(i,k,j) = vy(i,k,j) - dvdp(i,k,j) * &
( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth

end do
end do
end do
Expand All @@ -135,14 +158,26 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
ptx(i,k,j) = ( ptg(ip,k,j) - ptg(im,k,j) ) / (2._r8 * grid%dl)
ptx(i,k,j) = ptx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))

! Topography correction
ptx(i,k,j) = ptx(i,k,j) - dptdp(i,k,j) * &
( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))

! U-wind
ux(i,k,j) = ( ug(ip,k,j) - ug(im,k,j) ) / (2._r8 *grid%dl)
! Topography correction
ux(i,k,j) = ux(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))

! Topography correction
ux(i,k,j) = ux(i,k,j) - dudp(i,k,j) * &
( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))

Comment on lines +167 to +173
Copy link
Collaborator

@cacraigucar cacraigucar Dec 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wondering if you can give a bit more detail on "topography correction" or remove lines 169 and 170 (the blank line and duplicated comment) so the comment applies to both lines.

Copy link
Collaborator

@chihchen24 chihchen24 Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Horizontal temperature gradient for fronotogenesis function was originally computed over terrain-following coordinate. The correction terms are added so horizontal temperature gradient is on isobaric coordinate. This applies for all "correction" terms.

! V-wind
vx(i,k,j) = ( vg(ip,k,j) - vg(im,k,j) ) / (2._r8 *grid%dl)
vx(i,k,j) = vx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))

! Topography correction
vx(i,k,j) = vx(i,k,j) - dvdp(i,k,j) * &
( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))
end do
end do
end do
Expand Down Expand Up @@ -180,6 +215,66 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)

end subroutine gws_src_fnct

subroutine compute_vertical_derivative(grid,pint,pmid,data,ddata_dp)
!---------------------------------------------------------------------------
use dynamics_vars, only: t_fvdycore_grid

type (t_fvdycore_grid), intent(in) :: grid ! grid for XY decomp
real(r8), intent(in ) :: pint(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy) ! interface pressure
real(r8), intent(in ) :: pmid(grid%ifirstxy:grid%ilastxy ,plev,grid%jfirstxy:grid%jlastxy) ! mid-point pressure
real(r8), intent(in ) :: data(grid%ifirstxy:grid%ilastxy ,plev,grid%jfirstxy:grid%jlastxy)
real(r8), intent(out) :: ddata_dp(grid%ifirstxy:grid%ilastxy ,plev,grid%jfirstxy:grid%jlastxy)
!---------------------------------------------------------------------------
integer :: i,j,k
real(r8) :: pint_above ! pressure interpolated to interface above the current k mid-point
real(r8) :: pint_below ! pressure interpolated to interface below the current k mid-point
real(r8) :: dint_above ! data interpolated to interface above the current k mid-point
real(r8) :: dint_below ! data interpolated to interface below the current k mid-point
integer :: beglatxy, endlatxy, beglonxy, endlonxy
!---------------------------------------------------------------------------
beglatxy = grid%jfirstxy
endlatxy = grid%jlastxy
beglonxy = grid%ifirstxy
endlonxy = grid%ilastxy

!$omp parallel do private (i,j,k,pint_above,pint_below,dint_above,dint_below)
do j = beglatxy, endlatxy
do k = 2, plev-1
do i = beglonxy, endlonxy
pint_above = pint(i,k,j)
pint_below = pint(i,k+1,j)
dint_above = ( data(i,k-1,j) + data(i,k,j) ) / 2.0_r8
dint_below = ( data(i,k+1,j) + data(i,k,j) ) / 2.0_r8
ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
end do
end do
end do

k = 1
!$omp parallel do private (i,j,pint_above,pint_below,dint_above,dint_below)
do j = beglatxy, endlatxy
do i = beglonxy, endlonxy
pint_above = pmid(i,k,j)
pint_below = pint(i,k+1,j)
dint_above = data(i,k,j)
dint_below = ( data(i,k+1,j) + data(i,k,j) ) / 2.0_r8
ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
end do
end do

k = plev
!$omp parallel do private (i,j,pint_above,pint_below,dint_above,dint_below)
do j = beglatxy, endlatxy
do i = beglonxy, endlonxy
pint_above = pint(i,k,j)
pint_below = pmid(i,k,j)
dint_above = ( data(i,k-1,j) + data(i,k,j) ) / 2.0_r8
dint_below = data(i,k,j)
ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
end do
end do
end subroutine compute_vertical_derivative

subroutine ghost_array(grid, x, xg)

! subroutine added by Francis Vitt -- 7 July 2008
Expand Down
8 changes: 2 additions & 6 deletions src/dynamics/se/dycore/global_norms_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -578,7 +578,7 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,&

call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_p ,1.0_r8 ,'_p ')
call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu ,1.0_r8,' ')
call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,2.5_r8 ,'_div')
call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,1.0_r8 ,'_div')

if (nu_q<0) nu_q = nu_p ! necessary for consistency
if (nu_t<0) nu_t = nu_p ! temperature damping is always equal to nu_p
Expand Down Expand Up @@ -639,11 +639,7 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,&
!
! if user or namelist is not specifying sponge del4 settings here are best guesses (empirically determined)
!
if (top_042_090km) then
if (sponge_del4_lev <0) sponge_del4_lev = 4
if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 3.375_r8 !max value without having to increase subcycling of div4
if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 3.375_r8 !max value without having to increase subcycling of div4
else if (top_090_140km.or.top_140_600km) then ! defaults for waccm(x)
if (top_090_140km.or.top_140_600km) then ! defaults for waccm(x)
if (sponge_del4_lev <0) sponge_del4_lev = 20
if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 5.0_r8
if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 10.0_r8
Expand Down
13 changes: 8 additions & 5 deletions src/dynamics/se/dycore/prim_state_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
use se_dyn_time_mod, only: tstep
use control_mod, only: rsplit, qsplit
use perf_mod, only: t_startf, t_stopf
use hycoef, only: hyai, ps0
type (element_t), intent(inout) :: elem(:)
type (TimeLevel_t), target, intent(in) :: tl
type (hybrid_t), intent(in) :: hybrid
Expand Down Expand Up @@ -62,7 +63,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
! moist surface pressure
if (use_cslam) then
do ie=nets,nete
moist_ps_fvm(:,:,ie)=SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
moist_ps_fvm(:,:,ie)=SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)+hyai(1)*ps0
do q=dry_air_species_num+1,thermodynamic_active_species_num
m_cnst = thermodynamic_active_species_idx(q)
do k=1,nlev
Expand Down Expand Up @@ -134,8 +135,8 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
max_local(ie,5) = 0.0_r8
end if
if (use_cslam) then
min_local(ie,6) = MINVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))
max_local(ie,6) = MAXVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))
min_local(ie,6) = MINVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))+hyai(1)*ps0
max_local(ie,6) = MAXVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))+hyai(1)*ps0
min_local(ie,7) = MINVAL(moist_ps_fvm(:,:,ie))
max_local(ie,7) = MINVAL(moist_ps_fvm(:,:,ie))
min_local(ie,8) = MINVAL(elem(ie)%state%psdry(:,:))
Expand Down Expand Up @@ -207,7 +208,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%c(1:nc,1:nc,:,q)*fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
end do
q=statediag_numtrac+1
tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)+hyai(1)*ps0
q=statediag_numtrac+2
tmp_fvm(:,:,q,ie) = moist_ps_fvm(:,:,ie)
end do
Expand Down Expand Up @@ -291,7 +292,9 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
101 format (A12,A23,A23,A23,A23)

#ifdef waccm_debug
call prim_printstate_cslam_gamma(elem, tl,hybrid,nets,nete, fvm)
if (use_cslam) then
call prim_printstate_cslam_gamma(elem, tl,hybrid,nets,nete, fvm)
end if
#endif
call prim_printstate_U(elem, tl,hybrid,nets,nete, fvm)
end subroutine prim_printstate
Expand Down