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

Jump to bottom

Phase 2 of GW development #1117

Draft
wants to merge 2 commits into
base: cam_development
Choose a base branch
from
Draft

Phase 2 of GW development #1117

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dd84fa1
phase 2 of gw devel
JulioTBacmeister Aug 5, 2024
176bbd3
add template for two namelist variables for gw_drag
cacraigucar Aug 23, 2024
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
5 changes: 5 additions & 0 deletions bld/build-namelist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3795,6 +3795,7 @@ my $do_gw_convect_sh = ($nl->get_value('use_gw_convect_sh') =~ /$TRUE/io);
my $do_gw_movmtn_pbl = ($nl->get_value('use_gw_movmtn_pbl') =~ /$TRUE/io);
my $do_gw_rdg_beta = ($nl->get_value('use_gw_rdg_beta') =~ /$TRUE/io);
my $do_gw_rdg_gamma = ($nl->get_value('use_gw_rdg_gamma') =~ /$TRUE/io);
my $do_gw_rdg_resid = ($nl->get_value('use_gw_rdg_resid') =~ /$TRUE/io);

my $do_divstream = ($nl->get_value('gw_rdg_do_divstream') =~ /$TRUE/io);

Expand Down Expand Up @@ -3872,6 +3873,10 @@ if ($do_gw_rdg_beta) {
add_default($nl, 'gw_prndl');
}

if ($do_gw_rdg_resid) {
add_default($nl, 'effgw_rdg_resid');
}

if ($do_gw_rdg_gamma) {
add_default($nl, 'n_rdg_gamma', 'val'=>'-1');
add_default($nl, 'effgw_rdg_gamma', 'val'=>'1.0D0');
Expand Down
4 changes: 4 additions & 0 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 @@ -836,6 +836,10 @@
<effgw_rdg_beta_max model_top="lt" >0.5D0</effgw_rdg_beta_max>
<effgw_rdg_beta_max model_top="mt" >0.5D0</effgw_rdg_beta_max>

<effgw_rdg_beta_max >JULIO</effgw_rdg_beta_max>
<effgw_rdg_beta_max model_top="lt" >JULIO</effgw_rdg_beta_max>
<effgw_rdg_beta_max model_top="mt" >JULIO</effgw_rdg_beta_max>

<!-- setting for gravity waves from shallow convection. -->
<effgw_beres_sh>0.03D0</effgw_beres_sh>

Expand Down
12 changes: 12 additions & 0 deletions bld/namelist_files/namelist_definition.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1338,6 +1338,12 @@ Whether or not to enable gravity waves from PBL moving mountains source.
Default: .false.
</entry>

<entry id="use_gw_rdg_resid" type="logical" category="gw_drag"
group="phys_ctl_nl" valid_values="" >
JULIO - need to add comment
Default: .false.
</entry>

<entry id="pgwv" type="integer" category="gw_drag"
group="gw_drag_nl" valid_values="" >
Gravity wave spectrum dimension (wave numbers are from -pgwv to pgwv).
Expand Down Expand Up @@ -1441,6 +1447,12 @@ Max efficiency associated with anisotropic OGW.
Default: 1.0
</entry>

<entry id="effgw_rdg_resid" type="real" category="gw_drag"
group="gw_drag_nl" valid_values="" >
JULIO - need to add comment and appropriate default
Default: ?????? JULIO
</entry>

<entry id="rdg_beta_cd_llb" type="real" category="gw_drag"
group="gw_drag_nl" valid_values="" >
Drag coefficient for obstacles in low-level flow.
Expand Down
124 changes: 119 additions & 5 deletions src/physics/cam/gw_drag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,11 @@ module gw_drag
! Beres (shallow convection).
real(r8) :: effgw_beres_sh = unset_r8


! JULIO - Please put in appropriate comment
logical :: use_gw_rdg_resid = false
read(r8) :: effgw_rdg_resid = unset_r8

! Horzontal wavelengths [m].
real(r8), parameter :: wavelength_mid = 1.e5_r8
real(r8), parameter :: wavelength_long = 3.e5_r8
Expand Down Expand Up @@ -168,7 +173,9 @@ module gw_drag
integer, parameter :: prdg = 16

real(r8), allocatable, dimension(:,:), target :: &
rdg_gbxar
rdg_gbxar, &
rdg_isovar, &
rdg_isowgt

! Meso Beta
real(r8), allocatable, dimension(:,:,:), target :: &
Expand Down Expand Up @@ -250,7 +257,9 @@ subroutine gw_drag_readnl(nlfile)
rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, bnd_rdggm, &
gw_oro_south_fac, gw_limit_tau_without_eff, &
gw_lndscl_sgh, gw_prndl, gw_apply_tndmax, gw_qbo_hdepth_scaling, &
gw_top_taper, front_gaussian_width, alpha_gw_movmtn
gw_top_taper, front_gaussian_width, alpha_gw_movmtn, use_gw_rdg_resid, &
effgw_rdg_resid

!----------------------------------------------------------------------

if (use_simple_phys) return
Expand Down Expand Up @@ -359,6 +368,11 @@ subroutine gw_drag_readnl(nlfile)
call mpi_bcast(alpha_gw_movmtn, 1, mpi_real8, mstrid, mpicom, ierr)
if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: alpha_gw_movmtn")

call mpi_bcast(use_gw_rdg_resid, 1, mpi_logical, mstrid, mpicom, ierr)
if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_resid")
call mpi_bcast(effgw_rdg_resid, 1, mpi_real8, mstrid, mpicom, ierr)
if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_resid")

! Check if fcrit2 was set.
call shr_assert(fcrit2 /= unset_r8, &
"gw_drag_readnl: fcrit2 must be set via the namelist."// &
Expand Down Expand Up @@ -648,6 +662,8 @@ subroutine gw_init()
! Get beta ridge data
allocate( &
rdg_gbxar(pcols,begchunk:endchunk), &
rdg_isovar(pcols,begchunk:endchunk), &
rdg_isowgt(pcols,begchunk:endchunk), &
rdg_hwdth(pcols,prdg,begchunk:endchunk), &
rdg_clngt(pcols,prdg,begchunk:endchunk), &
rdg_mxdis(pcols,prdg,begchunk:endchunk), &
Expand All @@ -659,6 +675,14 @@ subroutine gw_init()
if (.not. found) call endrun(sub//': ERROR: GBXAR not found on topo file')
rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2

call infld('ISOVAR', fh_topo, dim1name, dim2name, 1, pcols, &
begchunk, endchunk, rdg_isovar, found, gridname='physgrid')
if (.not. found) call endrun(sub//': ERROR: ISOVAR not found on topo file')

call infld('ISOWGT', fh_topo, dim1name, dim2name, 1, pcols, &
begchunk, endchunk, rdg_isowgt, found, gridname='physgrid')
if (.not. found) call endrun(sub//': ERROR: ISOWGT not found on topo file')

call infld('HWDTH', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
1, prdg, begchunk, endchunk, rdg_hwdth, found, gridname='physgrid')
if (.not. found) call endrun(sub//': ERROR: HWDTH not found on topo file')
Expand Down Expand Up @@ -734,15 +758,39 @@ subroutine gw_init()
call addfld('ZMGW', (/ 'lev' /) , 'A' ,'m' , &
'midlayer geopotential heights in GW code ' )


call addfld('NIEGW', (/ 'ilev' /) , 'I' ,'1/s' , &
'interface BV freq in GW code ' )
call addfld('NMEGW', (/ 'lev' /) , 'I' ,'1/s' , &
'midlayer BV freq in GW code ' )
call addfld('RHOIEGW', (/ 'ilev' /) , 'I' ,'kg/m^3' , &
'interface density in GW code ' )
call addfld('PINTEGW', (/ 'ilev' /) , 'I' ,'Pa' , &
'interface density in GW code ' )

call addfld('TAUM1_DIAG' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')
call addfld('TAU1RDGBETAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')
call addfld('UBM1BETA', (/ 'lev' /) , 'A' ,'s-1' , &
call addfld('UBM1BETA', (/ 'lev' /) , 'A' ,'m s-1' , &
'On-ridge wind profile ' )
call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I' ,'m s-1' , &
call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I' ,'m s-2' , &
'On-ridge wind tendency from ridge 1 ')

call addfld('TAURESIDBETAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')
call addfld('UBMRESIDBETA', (/ 'lev' /) , 'I' ,'m s-1' , &
'On-ridge wind profile ' )
call addfld('UBIRESIDBETA', (/ 'ilev' /) , 'I' ,'m s-1' , &
'On-ridge wind profile (interface) ' )
call addfld('SRC_LEVEL_RESIDBETA', horiz_only , 'I' ,'1' , &
'src level index for ridge residual ' )
call addfld('TAUORO_RESID', horiz_only , 'I' ,'N m-2' , &
'Surface mom flux from ridge reisdual ' )
call addfld('TAUDIAG_RESID' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')


do i = 1, 6
write(cn, '(i1)') i
call addfld('TAU'//cn//'RDGBETAY' , (/ 'ilev' /), 'I', 'N m-2', &
Expand Down Expand Up @@ -1580,6 +1628,12 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
real(r8), pointer :: angll(:,:)
! anisotropy of ridges.
real(r8), pointer :: anixy(:,:)
! sqrt(residual variance) not repr by ridges (assumed isotropic).
real(r8), pointer :: isovar(:)
! area fraction of res variance
real(r8), pointer :: isowgt(:)



! Gamma ridges
! width of ridges.
Expand Down Expand Up @@ -2257,6 +2311,8 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
mxdis => rdg_mxdis(:ncol,:,lchnk)
angll => rdg_angll(:ncol,:,lchnk)
anixy => rdg_anixy(:ncol,:,lchnk)
isovar => rdg_isovar(:ncol,lchnk)
isowgt => rdg_isowgt(:ncol,lchnk)

where(mxdis < 0._r8)
mxdis = 0._r8
Expand All @@ -2276,6 +2332,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
nm, ni, rhoi, kvtt, q, dse, &
effgw_rdg_beta, effgw_rdg_beta_max, &
hwdth, clngt, gbxar, mxdis, angll, anixy, &
isovar, isowgt, &
rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
ptend, flx_heat)

Expand Down Expand Up @@ -2306,6 +2363,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
nm, ni, rhoi, kvtt, q, dse, &
effgw_rdg_gamma, effgw_rdg_gamma_max, &
hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, &
isovar, isowgt, &
rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, &
ptend, flx_heat)

Expand Down Expand Up @@ -2347,11 +2405,12 @@ subroutine gw_rdg_calc( &
effgw_rdg, effgw_rdg_max, &
hwdth, clngt, gbxar, &
mxdis, angll, anixy, &
isovar, isowgt, &
rdg_cd_llb, trpd_leewv, &
ptend, flx_heat)

use coords_1d, only: Coords1D
use gw_rdg, only: gw_rdg_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
use gw_rdg, only: gw_rdg_src, gw_rdg_resid_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
use gw_common, only: gw_drag_prof, energy_change

character(len=5), intent(in) :: type ! BETA or GAMMA
Expand Down Expand Up @@ -2385,6 +2444,9 @@ subroutine gw_rdg_calc( &
real(r8), intent(in) :: angll(ncol,prdg) ! orientation of ridges.
real(r8), intent(in) :: anixy(ncol,prdg) ! Anisotropy parameter.

real(r8), intent(in) :: isovar(ncol) ! sqrt of residual variance
real(r8), intent(in) :: isowgt(ncol) ! area frac of residual variance

real(r8), intent(in) :: rdg_cd_llb ! Drag coefficient for low-level flow
logical, intent(in) :: trpd_leewv

Expand Down Expand Up @@ -2604,6 +2666,58 @@ subroutine gw_rdg_calc( &

end do ! end of loop over multiple ridges

! Add additional GW from residual variance. Assumed isotropic
!kwvrdg = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!!
kwvrdg = 0.001_r8 / ( 100._r8 )
effgw = 1.0_r8 * isowgt
tauoro = 0._r8

call gw_rdg_resid_src(ncol, band_oro, p, &
u, v, t, isovar, kwvrdg, zi, nm, &
src_level, tend_level, tau, ubm, ubi, xv, yv, &
ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, phase_speeds, tauoro )

call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
t, vramp, &
piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
kwvrdg=kwvrdg, &
satfac_in = 1._r8, lapply_vdiff=gw_rdg_do_vdiff , tau_diag=tau_diag )

! Add the tendencies from isotropic residual to the totals.
do k = 1, pver
! diagnostics
utrdg(:,k) = utrdg(:,k) + utgw(:,k)
vtrdg(:,k) = vtrdg(:,k) + vtgw(:,k)
ttrdg(:,k) = ttrdg(:,k) + ttgw(:,k)
! physics tendencies
ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
end do

do m = 1, pcnst
do k = 1, pver
ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
end do
end do

do k = 1, pver+1
taurx0(:,k) = tau(:,0,k)*xv
taury0(:,k) = tau(:,0,k)*yv
taurx(:,k) = taurx(:,k) + taurx0(:,k)
taury(:,k) = taury(:,k) + taury0(:,k)
end do

call outfld('TAUDIAG_RESID', tau_diag, ncol, lchnk)
call outfld('TAUORO_RESID', tauoro , ncol, lchnk)
call outfld('TAURESID'//trim(type)//'M', tau(:,0,:), ncol, lchnk)
call outfld('UBMRESID'//trim(type), ubm, ncol, lchnk)
call outfld('UBIRESID'//trim(type), ubi, ncol, lchnk)
call outfld('SRC_LEVEL_RESID'//trim(type), 1._r8*src_level , ncol, lchnk)
! end of residual variance calc

! Calculate energy change for output to CAM's energy checker.
call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
ptend%v(:ncol,:), ptend%s(:ncol,:), de)
Expand Down
Loading
Loading