Skip to content

Commit

Permalink
revert some fv dycore changes for ssfcst stability
Browse files Browse the repository at this point in the history
  • Loading branch information
@TeaganKing
TeaganKing committed Feb 9, 2024
1 parent a03b84b commit 9e3e419
Show file tree
Hide file tree
Showing 9 changed files with 80 additions and 148 deletions.
119 changes: 34 additions & 85 deletions src/dynamics/fv/cd_core.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
mlt, ncx, ncy, nmfx, nmfy, iremote, &
cxtag, cytag, mfxtag, mfytag, &
cxreqs, cyreqs, mfxreqs, mfyreqs, &
kmtp, am_correction, am_geom_crrct, am_fixer, &
dod, don, high_order_top)
kmtp, am_correction, am_fixer, dod, don ,high_order_top)

! Dynamical core for both C- and D-grid Lagrangian dynamics
!
Expand Down Expand Up @@ -83,7 +82,6 @@ subroutine cd_core(grid, nx, u, v, pt, &
real(r8), intent(in) :: del2coef
integer, intent(in) :: kmtp ! range of levels (1:kmtp) where order is reduced
logical, intent(in) :: am_correction ! logical switch for correction (applied here)
logical, intent(in) :: am_geom_crrct ! logical switch for correction (applied here)
logical, intent(in) :: am_fixer ! logical switch for fixer (generate out args)
logical, intent(in) :: high_order_top ! use uniform 4th order everywhere (incl. model top)

Expand Down Expand Up @@ -170,7 +168,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
real(r8), intent(out) :: &
ptk(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)

! Omega calculation
! C.-C. Chen, omega calculation
real(r8), intent(out) :: &
cx_om(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast) ! Courant in X
real(r8), intent(out) :: &
Expand Down Expand Up @@ -231,7 +229,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
integer :: npes_yz

integer i, j, k, ml
integer js1g1, js2g0, js2g1, jn2g1, js4g0, jn3g0
integer js1g1, js2g0, js2g1, jn2g1 ,js4g0,jn3g0
integer jn2g0, jn1g1
integer iord , jord

Expand Down Expand Up @@ -291,10 +289,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
! Used to allow the same code to execute with or without the AM correction
real(r8) :: ptr(grid%im,grid%jfirst-1:grid%jlast+1,grid%kfirst:grid%klast+1)

logical :: sw_am_corr
logical :: am_press_crrct
real(r8) :: wg_hiord
real(r8) :: tpr, acap
logical :: sw_am_corr

!******************************************************************
!******************************************************************
Expand Down Expand Up @@ -376,13 +371,6 @@ subroutine cd_core(grid, nx, u, v, pt, &
endif
#endif

am_press_crrct = am_geom_crrct .and. am_correction
if (am_press_crrct) then
wg_hiord = -1._r8
else
wg_hiord = 0._r8
endif

npes_yz = grid%npes_yz

im = grid%im
Expand Down Expand Up @@ -411,21 +399,16 @@ subroutine cd_core(grid, nx, u, v, pt, &
kelp(grid%im,grid%jfirst-1:grid%jlast ,grid%kfirst:grid%klast ), &
dpn(grid%im,grid%jfirst :grid%jlast ,grid%kfirst:grid%klast ), &
dpo(grid%im,grid%jfirst :grid%jlast ,grid%kfirst:grid%klast ) )
acap = 1._r8/4._r8 ! effective AM/MoI contribution from polar caps
endif
if (am_press_crrct) then
allocate( &
dpr(grid%im,grid%jfirst-1:grid%jlast+1,grid%kfirst:grid%klast ) )
xakap = 1._r8/cap3vc(1,jfirst,kfirst)
endif
if (am_correction) then
allocate( &
dpr(grid%im,grid%jfirst-1:grid%jlast+1,grid%kfirst:grid%klast ), &
ddpu(grid%im,grid%jfirst :grid%jlast ,grid%kfirst:grid%klast ), &
dpns(grid%jfirst:grid%jlast,grid%kfirst:grid%klast), &
ddus(grid%jfirst:grid%jlast,grid%kfirst:grid%klast) )
ddus = 0._r8
xakap = 1._r8/cap3vc(1,jfirst,kfirst)
else
xakap = 1._r8
xakap = 1._r8
endif

! maintain consistent accuracy (uniform PPM order) over domain
Expand Down Expand Up @@ -471,7 +454,6 @@ subroutine cd_core(grid, nx, u, v, pt, &
jn2g0 = min(jm-1,jlast)
jn1g1 = min(jm,jlast+1)
jn2g1 = min(jm-1,jlast+1)

js4g0 = max(4,jfirst)
jn3g0 = min(jm-2,jlast)

Expand Down Expand Up @@ -750,7 +732,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
ua(1,jfirst-ng_d,k), va(1,jfirst-ng_s,k), &
uc(1,jfirst-ng_d,k), vc(1,jfirst-2,k), &
u_cen(1,jfirst-ng_d,k), v_cen(1,jfirst-ng_s,k), &
reset_winds, met_rlx(k), am_geom_crrct)
reset_winds, met_rlx(k), am_correction)

! Optionally filter advecting C-grid winds
if (filtcw .gt. 0) then
Expand Down Expand Up @@ -798,7 +780,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
ua(1,jfirst-ng_d,k), va(1,jfirst-ng_s,k), &
uc(1,jfirst-ng_d,k), vc(1,jfirst-2,k), &
ptc(1,jfirst,k), delpf(1,jfirst-ng_d,k), &
ptk(1,jfirst,k), tiny, iord, jord, am_geom_crrct)
ptk(1,jfirst,k), tiny, iord, jord, am_correction)
end do

call FVstopclock(grid,'---C_CORE')
Expand Down Expand Up @@ -1062,7 +1044,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
do k = kfirst, klast
do j = js2g0, jn2g0

if (am_press_crrct) then
if (am_correction) then

do i = 1, im
! AM fix: ensure interior pressure torque vanishes
Expand Down Expand Up @@ -1143,7 +1125,7 @@ subroutine cd_core(grid, nx, u, v, pt, &

call FVstartclock(grid,'---C_V_PGRAD')

if (am_press_crrct) then
if (am_correction) then
!$omp parallel do private(i, j, k)
! AM correction (pressure, advective winds): pxc -> ptr
do k = kfirst, klast+1
Expand Down Expand Up @@ -1326,20 +1308,8 @@ subroutine cd_core(grid, nx, u, v, pt, &
dpo(i,j,k)=(kelp(i,j,k)*cosp(j)+kelp(i,j-1,k)*cosp(j-1))/(cosp(j)+cosp(j-1)) ! A->D
end do
end do
if (jfirst == 1) then
do i = 1, im
dpn(i, 2,k)=(help(i, 2 ,k)*cosp( 2 )+acap*help(i, 1,k)*cose( 2))/cosp( 2 )
dpo(i, 2,k)=(kelp(i, 2 ,k)*cosp( 2 )+acap*kelp(i, 1,k)*cose( 2))/cosp( 2 )
end do
endif
if (jlast == jm) then
do i = 1, im
dpn(i,jm,k)=(help(i,jm-1,k)*cosp(jm-1)+acap*help(i,jm,k)*cose(jm))/cosp(jm-1)
dpo(i,jm,k)=(kelp(i,jm-1,k)*cosp(jm-1)+acap*kelp(i,jm,k)*cose(jm))/cosp(jm-1)
end do
endif
end do

if (am_correction) then
!$omp parallel do private(i, j, k)
do k = kfirst, klast
Expand All @@ -1354,30 +1324,23 @@ subroutine cd_core(grid, nx, u, v, pt, &
do k = kfirst, klast
do j = js2g0, jlast
do i = 1, im
ddu(i,j,k) = ddu(i,j,k)*0.5_r8*(dpo(i,j,k) + dpn(i,j,k))
ddu(i,j,k)=ddu(i,j,k)* D0_5*(dpo(i,j,k)+dpn(i,j,k)*3._r8)*D0_5
end do
end do
end do

!$omp parallel do private(i, j, k)
do k = kfirst, klast
do j = js2g0, jlast
ddus(j,k) = ddu(1,j,k) &
+ (u(1,j,k) + uc(1,j,k)*0.5_r8)*ddpu(1,j,k) &
+ wg_hiord*vf(1,j,k)*(dpn(1,j,k) - dpo(1,j,k))*0.5_r8
ddus(j,k) = ddu(1,j,k) + (u(1,j,k) + uc(1,j,k)/D4_0)*ddpu(1,j,k) - &
vf(1,j,k)*(dpn(1,j,k) - dpo(1,j,k))*D0_5
dpns(j,k) = dpn(1,j,k)
do i = 2, im
ddus(j,k) = ddus(j,k) &
+ ddu(i,j,k) &
+ (u(i,j,k)+uc(i,j,k)*0.5_r8)*ddpu(i,j,k) &
+ wg_hiord*vf(i,j,k)*(dpn(i,j,k)-dpo(i,j,k))*0.5_r8
ddus(j,k) = ddus(j,k) + ddu(i,j,k) +(u(i,j,k)+uc(i,j,k)/D4_0)*ddpu(i,j,k) - &
vf(i,j,k)*(dpn(i,j,k)-dpo(i,j,k))*D0_5
dpns(j,k) = dpns(j,k) + dpn(i,j,k)
end do
ddus(j,k) = ddus(j,k)/dpns(j,k)
! taper beyond 72S/N
tpr = max(abs(-2.5_r8 + ((j-1)-0.5_r8)*(5._r8/(jm-1))),2._r8)
tpr = cos(pi*tpr)**2
ddus(j,k)=ddus(j,k)*tpr
end do
end do

Expand All @@ -1393,29 +1356,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
end if ! (am_correction)

if (am_fixer) then
if (.not. am_geom_crrct) then
!$omp parallel do private(i, j, k)
do k = kfirst, klast
do j = js2g0, jlast
do i = 1, im
dpn(i,j,k) = (help(i,j,k) + help(i,j-1,k) )/ 2._r8
dpo(i,j,k) = (kelp(i,j,k) + kelp(i,j-1,k) )/ 2._r8
end do
end do
if (jfirst == 1) then
do i = 1, im
dpn(i,2,k) = (help(i,2,k) + help(i,1,k) )/ 2._r8
dpo(i,2,k) = (kelp(i,2,k) + kelp(i,1,k) )/ 2._r8
end do
endif
if (jlast == jm) then
do i = 1, im
dpn(i,jm,k) = (help(i,jm-1,k) + help(i,jm,k) )/ 2._r8
dpo(i,jm,k) = (kelp(i,jm-1,k) + kelp(i,jm,k) )/ 2._r8
end do
endif
end do
endif

!$omp parallel do private(i, j, k)
do k = kfirst, klast
do j = js2g0, jlast
Expand All @@ -1426,7 +1367,15 @@ subroutine cd_core(grid, nx, u, v, pt, &
dod(j,k) = dod(j,k) + (cosp(j) + cosp(j-1))*cose(j)**2*dpn(i,j,k)
end do
end do

! north pole
if (jfirst == 1) then
do i = 1, im
dod(1,k) = dod(1,k) + grid%acap/(D0_5*im)*cose(1)**2*help(i,1,k)
end do
end if
end do

end if ! (am_fixer)

call FVstopclock(grid,'---dp4corr_COMM_2')
Expand Down Expand Up @@ -1734,8 +1683,8 @@ subroutine cd_core(grid, nx, u, v, pt, &
end if
call FVstopclock(grid,'---PRE_D_PGRAD_COMM_1')
#endif
if (am_press_crrct) then

if (am_correction) then
! AM correction (pressure, prognostic winds): pkc -> ptr
!$omp parallel do private(i, j, k)
do k = kfirst, klast+1
Expand All @@ -1756,7 +1705,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
end do
endif

if (am_press_crrct) then
if (am_correction) then
!$omp parallel do private(i, j, k)
! Beware k+1 references directly below (AAM)
do k = kfirst, klast
Expand Down Expand Up @@ -1809,7 +1758,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
cycle
end if

if (am_press_crrct) then
if (am_correction) then
do j=js2g1,jn2g0 ! wk3 needed S
wk3(1,j) = (wz(1,j,k)+wz(im,j,k)) * &
(ptr(1,j,k) - ptr(im,j,k))
Expand Down Expand Up @@ -1859,7 +1808,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
! N-S walls

do j=js2g0,jn1g1 ! wk1 needed N
if (am_press_crrct) then
if (am_correction) then
do i=1,im
wk1(i,j) = (wz(i,j,k) + wz(i,j-1,k))*(ptr(i,j,k) - ptr(i,j-1,k))
enddo
Expand Down Expand Up @@ -1897,7 +1846,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
enddo
enddo

if (am_press_crrct) then
if (am_correction) then

! use true pressure for wk1, then update it
do j = js1g1, jn1g1
Expand All @@ -1924,7 +1873,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
!$omp parallel do private(i, j, k, wk, wk1, wk2, wk3)
do k = kfirst, klast

if (am_press_crrct) then
if (am_correction) then
do j = js1g1, jn1g1
wk1(1,j) = dpr(1,j,k) + dpr(im,j,k)
do i = 2, im
Expand Down Expand Up @@ -1964,7 +1913,7 @@ subroutine cd_core(grid, nx, u, v, pt, &
* (wk2(im,j)-wk2(1,j)+wz3(im,j,k+1)-wz3(im,j,k))
end do

if (am_geom_crrct) then
if (am_correction) then
! apply cos-weighted avg'ing
do j = js2g0, jn2g0 ! Assumes wk2 ghosted on N
do i = 1, im
Expand Down
8 changes: 4 additions & 4 deletions src/dynamics/fv/d2a3dikj.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module d2a3dikj_mod
!
! !INTERFACE:

subroutine d2a3dikj(grid, am_geom_crrct, u, v, ua, va)
subroutine d2a3dikj(grid, am_correction, u, v, ua, va)

! !USES:

Expand All @@ -36,7 +36,7 @@ subroutine d2a3dikj(grid, am_geom_crrct, u, v, ua, va)
implicit none
! !INPUT PARAMETERS:
type (t_fvdycore_grid), intent(in) :: grid
logical, intent(in) :: am_geom_crrct
logical, intent(in) :: am_correction
real(r8), intent(in) :: u(grid%ifirstxy:grid%ilastxy, &
grid%jfirstxy:grid%jlastxy,grid%km) ! U-Wind
real(r8), intent(in) :: v(grid%ifirstxy:grid%ilastxy, &
Expand Down Expand Up @@ -127,7 +127,7 @@ subroutine d2a3dikj(grid, am_geom_crrct, u, v, ua, va)

if ( jlastxy .lt. jm ) then

if (am_geom_crrct) then
if (am_correction) then
!$omp parallel do private(i, k)
do k = 1, km
do i = ifirstxy, ilastxy
Expand All @@ -146,7 +146,7 @@ subroutine d2a3dikj(grid, am_geom_crrct, u, v, ua, va)
end if
#endif

if (am_geom_crrct) then
if (am_correction) then
!$omp parallel do private(i,j,k)
do k = 1, km
do j = jfirstxy, jlastxy-1
Expand Down
10 changes: 5 additions & 5 deletions src/dynamics/fv/dp_coupling.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ subroutine d_p_coupling(grid, phys_state, phys_tend, pbuf2d, dyn_out)
allocate (v3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))

if (iam .lt. grid%npes_xy) then
call d2a3dikj(grid, dyn_state%am_geom_crrct, u3sxy, v3sxy, u3, v3)
call d2a3dikj(grid, dyn_state%am_correction, u3sxy, v3sxy, u3, v3)
end if ! (iam .lt. grid%npes_xy)

call t_stopf ('d2a3dikj')
Expand Down Expand Up @@ -303,9 +303,9 @@ subroutine d_p_coupling(grid, phys_state, phys_tend, pbuf2d, dyn_out)

if (iam .lt. grid%npes_xy) then
! (note dummy use of dva3 hence call order matters)
call d2a3dikj(grid, dyn_state%am_geom_crrct,duf3sxy, dummy, duf3 ,dva3)
call d2a3dikj(grid, dyn_state%am_geom_crrct,dua3sxy, dva3sxy, dua3, dva3)
call d2a3dikj(grid, dyn_state%am_geom_crrct, du3sxy, dv3sxy, du3 , dv3 )
call d2a3dikj(grid, dyn_state%am_correction,duf3sxy, dummy, duf3 ,dva3)
call d2a3dikj(grid, dyn_state%am_correction,dua3sxy, dva3sxy, dua3, dva3)
call d2a3dikj(grid, dyn_state%am_correction, du3sxy, dv3sxy, du3 , dv3 )
end if ! (iam .lt. grid%npes_xy)

call t_startf('DP_CPLN_fv_am')
Expand Down Expand Up @@ -950,7 +950,7 @@ subroutine p_d_coupling(grid, phys_state, phys_tend, &
call t_startf('uv3s_update')
if (iam .lt. grid%npes_xy) then
call uv3s_update(grid, dudtxy, u3sxy, dvdtxy, v3sxy, dt5, &
dyn_state%am_geom_crrct)
dyn_state%am_correction)
end if ! (iam .lt. grid%npes_xy)
call t_stopf('uv3s_update')

Expand Down
Loading

0 comments on commit 9e3e419

Please sign in to comment.