Skip to content
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

Add optional scaling to RRTMGP flux adjustment #668

Merged
merged 14 commits into from
Jul 9, 2021
Merged

Add optional scaling to RRTMGP flux adjustment #668

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
df30cf7
Limit full LW flux profile adjustment to below 100hPa.
dustinswales May 9, 2021
6961b10
Added more safeguards against out-of-bounds temperature to GP inputs.
dustinswales May 9, 2021
6ebe85e
Apply constant heating-rate adjustment above 100hPa
dustinswales May 11, 2021
230d479
Add vetical decay to impact of flux adjustment above threshold.
dustinswales May 11, 2021
5f7d697
Added logistic function to damp the LW flux adjustment with height
dustinswales May 24, 2021
2f6e708
Reorganized RRTMGP aerosol optics.
dustinswales May 24, 2021
7bc877d
Revert "Reorganized RRTMGP aerosol optics."
dustinswales May 25, 2021
eb7837d
Bug fix. Add transition to HR adjustment.
dustinswales May 25, 2021
f2d5570
Removed exp(1) from scaling
dustinswales May 26, 2021
ea01139
Housekeeping. Added comments.
dustinswales May 26, 2021
d93adbe
Further cleanup of dcyc2
dustinswales May 26, 2021
fe44d62
Revert "Added more safeguards against out-of-bounds temperature to GP…
dustinswales May 26, 2021
6a0e904
Omission from previous revert.
dustinswales May 27, 2021
41782f1
Change from PR review.
dustinswales May 27, 2021
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
86 changes: 56 additions & 30 deletions physics/dcyc2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ subroutine dcyc2t3_run &
& sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, &
& sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, &
& im, levs, deltim, fhswr, &
& dry, icy, wet, &
& dry, icy, wet, damp_LW_fluxadj, lfnc_k, lfnc_p0, &
& minGPpres, use_LW_jacobian, sfculw, fluxlwUP_jac, &
& t_lay, t_lev, p_lay, p_lev, flux2D_lwUP, flux2D_lwDOWN, &
& pert_radtend, do_sppt,ca_global, &
Expand Down Expand Up @@ -213,10 +213,11 @@ subroutine dcyc2t3_run &
! integer, intent(in) :: ipr
! logical lprnt
logical, dimension(:), intent(in) :: dry, icy, wet
logical, intent(in) :: use_LW_jacobian, pert_radtend
logical, intent(in) :: use_LW_jacobian, damp_LW_fluxadj, &
& pert_radtend
logical, intent(in) :: do_sppt,ca_global
real(kind=kind_phys), intent(in) :: solhr, slag, cdec, sdec, &
& deltim, fhswr, minGPpres
& deltim, fhswr, minGPpres, lfnc_k, lfnc_p0

real(kind=kind_phys), dimension(:), intent(in) :: &
& sinlat, coslat, xlon, coszen, tf, tsflw, sfcdlw, &
Expand Down Expand Up @@ -253,11 +254,19 @@ subroutine dcyc2t3_run &
integer, intent(out) :: errflg

! --- locals:
integer :: i, k, nstp, nstl, it, istsun(im),iSFC
integer :: i, k, nstp, nstl, it, istsun(im),iSFC,iTOA
real(kind=kind_phys) :: cns, coszn, tem1, tem2, anginc, &
& rstl, solang, dT
real(kind=kind_phys), dimension(im,levs+1) :: flxlwup_adj, &
& flxlwdn_adj, t_lev2
real(kind=kind_phys) :: fluxlwnet_adj,fluxlwnet,dT_sfc, &
&fluxlwDOWN_jac,lfnc,c1
! Length scale for flux-adjustment scaling
real(kind=kind_phys), parameter :: &
& L = 1.
! Scaling factor for downwelling LW Jacobian profile.
real(kind=kind_phys), parameter :: &
& gamma = 0.2
!
!===> ... begin here
!
Expand All @@ -267,9 +276,11 @@ subroutine dcyc2t3_run &

! Vertical ordering?
if (p_lev(1,1) .lt. p_lev(1, levs)) then
iSFC = levs
iSFC = levs + 1
Copy link
Collaborator

Choose a reason for hiding this comment

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

Was that an actual bug? Is iSFC for interfaces (which is always levs+1)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yes^2.

iTOA = 1
else
iSFC = 1
iTOA = levs + 1
endif

tem1 = fhswr / deltim
Expand Down Expand Up @@ -375,32 +386,47 @@ subroutine dcyc2t3_run &
call cmp_tlev(im, levs, minGPpres, p_lay, t_lay, p_lev, tsfc, &
& t_lev2)

! Compute adjusted net LW flux foillowing Hogan and Bozzo 2015 (10.1002/2015MS000455)
! Here we assume that the profile of the downwelling LW Jaconiam has the same shape
! as the upwelling, but scaled and offset.
! The scaling factor is 0.2
! The profile of the downwelling Jacobian (J) is offset so that
! J_dn_sfc / J_up_sfc = scaling_factor
! J_dn_toa / J_up_sfc = 0
!
! Adjust up/downward fluxes (at layer interfaces).
!
do k = 1, levs+1
do i = 1, im
dT = t_lev2(i,k) - t_lev(i,k)
flxlwup_adj(i,k) = flux2D_lwUP(i,k) + &
& fluxlwUP_jac(i,k)*dT
enddo
enddo
!
! Compute new heating rate (within each layer).
!
do k = 1, levs
htrlw(1:im,k) = &
& (flxlwup_adj(1:im,k+1) - flxlwup_adj(1:im,k) - &
& flux2D_lwDOWN(1:im,k+1) + flux2D_lwDOWN(1:im,k)) * &
& con_g / (con_cp * (p_lev(1:im,k+1) - p_lev(1:im,k)))
enddo

!
! Add radiative heating rates to physics heating rate
!
do k = 1, levs
do i = 1, im
dtdt(i,k) = dtdt(i,k) + swh(i,k)*xmu(i) + htrlw(i,k)
! Optionally, the flux adjustment can be damped with height using a logistic function
! fx ~ L / (1 + exp(-k*dp)), where dp = p - p0
! L = 1, fix scale between 0-1. - Fixed
! k = 1 / pressure decay length (Pa) - Controlled by namelist
! p0 = Transition pressure (Pa) - Controlled by namelsit
do i = 1, im
c1 = fluxlwUP_jac(i,iTOA) / fluxlwUP_jac(i,iSFC)
dT_sfc = t_lev2(i,iSFC) - t_lev(i,iSFC)
do k = 1, levs
! LW net flux
fluxlwnet = (flux2D_lwUP(i, k+1) - flux2D_lwUP(i, k) - &
& flux2D_lwDOWN(i,k+1) + flux2D_lwDOWN(i,k))
! Downward LW Jacobian (Eq. 9)
fluxlwDOWN_jac = gamma * &
& (fluxlwUP_jac(i,k)/fluxlwUP_jac(i,iSFC) - c1) / &
& (1 - c1)
! Adjusted LW net flux(Eq. 10)
fluxlwnet_adj = fluxlwnet + dT_sfc* &
& (fluxlwUP_jac(i,k)/fluxlwUP_jac(i,iSFC) - &
& fluxlwDOWN_jac)
! Adjusted LW heating rate
htrlw(i,k) = fluxlwnet_adj * con_g / &
& (con_cp * (p_lev(i,k+1) - p_lev(i,k)))

! Add radiative heating rates to physics heating rate. Optionally, scaled w/ height
! using a logistic function
if (damp_LW_fluxadj) then
lfnc = L / (1+exp(-(p_lev(i,k) - lfnc_p0)/lfnc_k))
else
lfnc = 1.
endif
dtdt(i,k) = dtdt(i,k) + swh(i,k)*xmu(i) + &
& htrlw(i,k)*lfnc + (1.-lfnc)*hlw(i,k)
enddo
enddo
else
Expand Down
26 changes: 26 additions & 0 deletions physics/dcyc2.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -370,6 +370,32 @@
type = logical
intent = in
optional = F
[damp_LW_fluxadj]
standard_name = flag_to_damp_RRTMGP_LW_jacobian_flux_adjustment
long_name = logical flag to control RRTMGP LW calculation
units = flag
dimensions = ()
type = logical
intent = in
optional = F
[lfnc_k]
standard_name = transition_pressure_length_scale_for_flux_damping
long_name = depth of transition layer in logistic function for LW flux adjustment damping
units = Pa
dimensions = ()
type = real
kind = kind_phys
intent = in
optional = F
[lfnc_p0]
standard_name = transition_pressure_for_flux_damping
long_name = transition pressure for LW flux adjustment damping
units = Pa
dimensions = ()
type = real
kind = kind_phys
intent = in
optional = F
[sfculw]
standard_name = surface_upwelling_longwave_flux_on_radiation_time_step
long_name = total sky sfc upward lw flux
Expand Down