Skip to content

Commit

Permalink
refactor gw_convect to include stealth bug fix and expose namelist pa…
Browse files Browse the repository at this point in the history 
…rameters
  • Loading branch information
@whannah1
whannah1 committed Nov 15, 2024
1 parent c551ee0 commit dbc44a1
Show file tree
Hide file tree
Showing 2 changed files with 139 additions and 86 deletions.
189 changes: 120 additions & 69 deletions components/eam/src/physics/cam/gw_convect.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,9 @@ module gw_convect
! This module handles gravity waves from convection, and was extracted from
! gw_drag in May 2013.
!

use gw_utils, only: r8

use gw_common, only: pver, pgwv
use cam_logfile, only: iulog
use gw_utils, only: r8
use gw_common, only: pver, pgwv

implicit none
private
Expand All @@ -21,8 +20,8 @@ module gw_convect
! Dimension for mean wind in heating.
integer :: maxuh

! Index for level at 700 mb.
integer :: k700
! Index for level for storm/steering flow (usually 700 mb)
integer :: k_src_wind

! Table of source spectra.
real(r8), allocatable :: mfcc(:,:,:)
Expand All @@ -31,19 +30,20 @@ module gw_convect

!==========================================================================

subroutine gw_convect_init(k700_in, mfcc_in, errstring)
! Index at 700 mb.
integer, intent(in) :: k700_in
! Source spectra to keep as table.
real(r8), intent(in) :: mfcc_in(:,:,:)
! Report any errors from this routine.
character(len=*), intent(out) :: errstring

subroutine gw_convect_init( plev_src_wind, mfcc_in, errstring)
use ref_pres, only: pref_edge
real(r8), intent(in) :: plev_src_wind ! previously hardcoded to 70000._r8
real(r8), intent(in) :: mfcc_in(:,:,:) ! Source spectra to keep as table
character(len=*), intent(out) :: errstring ! Report any errors from this routine
integer :: ierr

errstring = ""

k700 = k700_in
do k = 0, pver
if ( pref_edge(k+1) < plev_src_wind ) k_src_wind = k+1
end do

if (masterproc) write (iulog,*) 'gw_convect: steering flow level = ',k_src_wind

! First dimension is maxh.
maxh = size(mfcc_in,1)
Expand All @@ -60,7 +60,9 @@ end subroutine gw_convect_init

subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
zm, src_level, tend_level, tau, ubm, ubi, xv, yv, c, &
hdepth, maxq0, CF, hdepth_scaling_factor)
hdepth, maxq0_out, maxq0_conversion_factor, hdepth_scaling_factor, &
hdepth_min, storm_speed_min, &
use_gw_convect_old)
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
!
Expand Down Expand Up @@ -90,11 +92,20 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
real(r8), intent(in) :: zm(ncol,pver)

! Heating conversion factor
real(r8), intent(in) :: CF
real(r8), intent(in) :: maxq0_conversion_factor

! Scaling factor for the heating depth
real(r8), intent(in) :: hdepth_scaling_factor

! minimum hdepth for for spectrum lookup table
real(r8), intent(in) :: hdepth_min

! minimum convective storm speed
real(r8), intent(in) :: storm_speed_min

! switch for restoring legacy method
logical, intent(in) :: use_gw_convect_old

! Indices of top gravity wave source level and lowest level where wind
! tendencies are allowed.
integer, intent(out) :: src_level(ncol)
Expand All @@ -110,19 +121,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
real(r8), intent(out) :: c(ncol,-pgwv:pgwv)

! Heating depth and maximum heating in each column.
real(r8), intent(out) :: hdepth(ncol), maxq0(ncol)
real(r8), intent(out) :: hdepth(ncol), maxq0_out(ncol)

!---------------------------Local Storage-------------------------------
! Column and level indices.
integer :: i, k

! Zonal/meridional wind at 700mb.
real(r8) :: u700(ncol), v700(ncol)
! Zonal/meridional source wind
real(r8) :: u_src(ncol), v_src(ncol)
! 3.14...
real(r8), parameter :: pi = 4._r8*atan(1._r8)

! Maximum heating rate.
real(r8) :: q0(ncol)
real(r8) :: maxq0(ncol)

! Bottom/top heating range index.
integer :: mini(ncol), maxi(ncol)
Expand All @@ -135,36 +146,36 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
! Source level tau for a column.
real(r8) :: tau0(-PGWV:PGWV)
! Speed of convective cells relative to storm.
integer :: CS(ncol)
integer :: storm_speed(ncol)
! Index to shift spectra relative to ground.
integer :: shift

! Heating rate conversion factor. Change to take the value from caller and controllable by namelist (to tune QBO)
! real(r8), parameter :: CF = 20._r8
! Averaging length.
real(r8), parameter :: AL = 1.0e5_r8
! fixed parameters (we may want to expose these in the namelist for tuning)
real(r8), parameter :: tau_avg_length = 1.0e5_r8 ! spectrum averaging length
real(r8), parameter :: heating_altitude_max = 20e3 ! max altitude to check heating (probably don't need this)

integer :: ndepth_pos
integer :: ndepth_tot

!----------------------------------------------------------------------
! Initialize tau array
!----------------------------------------------------------------------

tau = 0.0_r8
tau = 0.0_r8
hdepth = 0.0_r8
q0 = 0.0_r8
tau0 = 0.0_r8
maxq0 = 0.0_r8
tau0 = 0.0_r8

!------------------------------------------------------------------------
! Determine 700 mb layer wind and unit vectors, then project winds.
! Determine source layer wind and unit vectors, then project winds.
!------------------------------------------------------------------------

! Just use the 700 mb interface values for the source wind speed and
! direction (unit vector).

u700 = u(:,k700)
v700 = v(:,k700)
! source wind speed and direction
u_src = u(:,k_src_wind)
v_src = v(:,k_src_wind)

! Get the unit vector components and magnitude at the surface.
call get_unit_vector(u700, v700, xv, yv, ubi(:,k700))
call get_unit_vector(u_src, v_src, xv, yv, ubi(:,k_src_wind))

! Project the local wind at midpoints onto the source wind.
do k = 1, pver
Expand All @@ -184,33 +195,62 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
! which heating rate is continuously positive.
!-----------------------------------------------------------------------

! First find the indices for the top and bottom of the heating range.
! Find indices for the top and bottom of the heating range.
mini = 0
maxi = 0
do k = pver, 1, -1
do i = 1, ncol

if (use_gw_convect_old) then
!---------------------------------------------------------------------
! original version used in CAM4/5/6 and EAMv1/2/3
do k = pver, 1, -1
do i = 1, ncol
if (mini(i) == 0) then
! Detect if we are outside the maximum range (where z = 20 km).
if (zm(i,k) >= 20000._r8) then
mini(i) = k
maxi(i) = k
else
! First spot where heating rate is positive.
if (netdt(i,k) > 0.0_r8) mini(i) = k
end if
! Detect if we are outside the maximum range (where z = 20 km).
if (zm(i,k) >= heating_altitude_max) then
mini(i) = k
maxi(i) = k
else
! First spot where heating rate is positive.
if (netdt(i,k) > 0.0_r8) mini(i) = k
end if
else if (maxi(i) == 0) then
! Detect if we are outside the maximum range (z = 20 km).
if (zm(i,k) >= 20000._r8) then
maxi(i) = k
else
! First spot where heating rate is no longer positive.
if (.not. (netdt(i,k) > 0.0_r8)) maxi(i) = k
end if
! Detect if we are outside the maximum range (z = 20 km).
if (zm(i,k) >= heating_altitude_max) then
maxi(i) = k
else
! First spot where heating rate is no longer positive.
if (.not. (netdt(i,k) > 0.0_r8)) maxi(i) = k
end if
end if
end do
! When all done, exit
if (all(maxi /= 0)) exit
end do
end do
! When all done, exit
if (all(maxi /= 0)) exit
end do
!---------------------------------------------------------------------
else
!---------------------------------------------------------------------
! cleaner version that addresses bug in original where heating max and
! depth were too low whenever heating <=0 occurred in the middle of
! the heating profile (ex. at the melting level)
do i = 1, ncol
do k = pver, 1, -1
if ( zm(i,k) < heating_altitude_max ) then
if ( netdt(i,k) > 0.0_r8 ) then
! Set mini as first spot where heating rate is positive
if ( mini(i)==0 ) mini(i) = k
! Set maxi to current level
maxi(i) = k
end if
else
! above the max check if indices were found
if ( mini(i)==0 ) mini(i) = k
if ( maxi(i)==0 ) maxi(i) = k
end if
end do
end do
!---------------------------------------------------------------------
end if


! Heating depth in km.
hdepth = [ ( (zm(i,maxi(i))-zm(i,mini(i)))/1000._r8, i = 1, ncol ) ]
Expand All @@ -223,19 +263,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
! Maximum heating rate.
do k = minval(maxi), maxval(mini)
where (k >= maxi .and. k <= mini)
q0 = max(q0, netdt(:ncol,k))
maxq0 = max(maxq0, netdt(:ncol,k))
end where
end do

!output max heating rate in K/day
maxq0 = q0*24._r8*3600._r8
maxq0_out = maxq0*24._r8*3600._r8

! Multipy by conversion factor
q0 = q0 * CF
maxq0 = maxq0 * maxq0_conversion_factor

! Taking ubm at 700 mb to be the storm speed, find the cell speed where
! the storm speed is > 10 m/s.
CS = int(sign(max(abs(ubm(:,k700))-10._r8, 0._r8), ubm(:,k700)))
! Taking ubm at assumed source level to be the storm speed,
! find the cell speed where the storm speed is > storm_speed_min
storm_speed = int(sign(max(abs(ubm(:,k_src_wind))-storm_speed_min, 0._r8), ubm(:,k_src_wind)))

uh = 0._r8
do k = minval(maxi), maxval(mini)
Expand All @@ -244,7 +284,7 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
end where
end do

uh = uh - real(CS, r8)
uh = uh - real(storm_speed, r8)

! Limit uh to table range.
uh = min(uh, real(maxuh, r8))
Expand All @@ -270,8 +310,19 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
!---------------------------------------------------------------------
! Look up spectrum only if depth >= 2.5 km, else set tau0 = 0.
!---------------------------------------------------------------------

if ((hdepth(i) >= 2.5_r8) .and. (abs(lat(i)) < (pi/2._r8))) then
if ((hdepth(i) >= hdepth_min)) then
ndepth_pos = 0
ndepth_tot = 0
do k = 1,pver
if ( k>=maxi(i).and.k<=mini(i) ) then
ndepth_tot = ndepth_tot + 1
if ( netdt(i,k)>0 ) ndepth_pos = ndepth_pos + 1
end if
end do
! write (iulog,*) 'WHDEBUG - i: ',i,' Hd: ',hdepth(i),' Hn: ',ndepth_pos,' N: ',ndepth_tot
end if

if ((hdepth(i) >= hdepth_min) .and. (abs(lat(i)) < (pi/2._r8))) then

!------------------------------------------------------------------
! Look up the spectrum using depth and uh.
Expand All @@ -280,11 +331,11 @@ subroutine gw_beres_src(ncol, ngwv, lat, u, v, netdt, &
tau0 = mfcc(NINT(hdepth(i)),NINT(uh(i)),:)

! Shift spectrum so that it is relative to the ground.
shift = -nint(real(CS(i), r8)/dc)
shift = -nint(storm_speed(i)/dc)
tau0 = cshift(tau0,shift)

! Adjust magnitude.
tau0 = tau0*q0(i)*q0(i)/AL
tau0 = tau0*maxq0(i)*maxq0(i)/tau_avg_length

! Adjust for critical level filtering.
Umini = max(nint(Umin(i)/dc),-PGWV)
Expand Down
36 changes: 19 additions & 17 deletions components/eam/src/physics/cam/gw_drag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,11 @@ module gw_drag
! namelist
logical :: history_amwg ! output the variables used by the AMWG diag package

logical :: use_gw_convect_old ! switch to enable legacy behavior
real(r8) :: gw_convect_plev_src_wind ! reference pressure level for source wind for convective GWD
real(r8) :: gw_convect_hdepth_min ! minimum hdepth for for convective GWD spectrum lookup table
real(r8) :: gw_convect_storm_speed_min ! minimum convective storm speed for convective GWD

!==========================================================================
contains
!==========================================================================
Expand All @@ -141,7 +146,9 @@ subroutine gw_drag_readnl(nlfile)

namelist /gw_drag_nl/ pgwv, gw_dc, tau_0_ubc, effgw_beres, effgw_cm, &
effgw_oro, fcrit2, frontgfc, gw_drag_file, taubgnd, gw_convect_hcf, &
hdepth_scaling_factor
hdepth_scaling_factor, gw_convect_hdepth_min, &
gw_convect_storm_speed_min, gw_convect_plev_src_wind, &
use_gw_convect_old)
!----------------------------------------------------------------------

if (masterproc) then
Expand Down Expand Up @@ -170,8 +177,12 @@ subroutine gw_drag_readnl(nlfile)
call mpibcast(frontgfc, 1, mpir8, 0, mpicom)
call mpibcast(taubgnd, 1, mpir8, 0, mpicom)
call mpibcast(gw_drag_file, len(gw_drag_file), mpichar, 0, mpicom)
call mpibcast(gw_convect_hcf, 1, mpir8, 0, mpicom)
call mpibcast(hdepth_scaling_factor, 1, mpir8, 0, mpicom)
call mpibcast(gw_convect_hcf, 1, mpir8, 0, mpicom)
call mpibcast(hdepth_scaling_factor, 1, mpir8, 0, mpicom)
call mpibcast(gw_convect_hdepth_min, 1, mpir8, 0, mpicom)
call mpibcast(gw_convect_storm_speed_min, 1, mpir8, 0, mpicom)
call mpibcast(gw_convect_plev_src_wind, 1, mpir8, 0, mpicom)
call mpibcast(use_gw_convect_old, 1, mpilog, 0, mpicom)
#endif

dc = gw_dc
Expand Down Expand Up @@ -224,7 +235,6 @@ subroutine gw_init()

! Index for levels at specific pressures.
integer :: kfront
integer :: k700

! output tendencies and state variables for CAM4 temperature,
! water vapor, cloud ice and cloud liquid budgets.
Expand Down Expand Up @@ -451,20 +461,10 @@ subroutine gw_init()

ttend_dp_idx = pbuf_get_index('TTEND_DP')

do k = 0, pver
! 700 hPa index
if (pref_edge(k+1) < 70000._r8) k700 = k+1
end do

if (masterproc) then
write (iulog,*) 'K700 =',k700
end if

! Initialization of Beres' parameterization parameters
call gw_init_beres(mfcc)
call gw_convect_init(k700, mfcc, errstring)
if (trim(errstring) /= "") &
call endrun("gw_convect_init: "//errstring)
call gw_convect_init(gw_convect_plev_src_wind, mfcc, errstring)
if (trim(errstring) /= "") call endrun("gw_convect_init: "//errstring)

! Output for gravity waves from the Beres scheme.
call gw_spec_addflds(prefix=beres_pf, scheme="Beres", &
Expand Down Expand Up @@ -765,7 +765,9 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)
! Determine wave sources for Beres04 scheme
call gw_beres_src(ncol, pgwv, state1%lat(:ncol), u, v, ttend_dp, &
zm, src_level, tend_level, tau, ubm, ubi, xv, yv, c, &
hdepth, maxq0, gw_convect_hcf, hdepth_scaling_factor)
hdepth, maxq0, gw_convect_hcf, hdepth_scaling_factor, &
gw_convect_hdepth_min, gw_convect_storm_speed_min, &
use_gw_convect_old)

do_latitude_taper = .false.

Expand Down

0 comments on commit dbc44a1

Please sign in to comment.