From dd84fa1b560bef88aa42635025b764a01958418d Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Mon, 5 Aug 2024 06:55:19 -0600 Subject: [PATCH 1/6] phase 2 of gw devel --- src/physics/cam/gw_drag.F90 | 110 ++++++++++++++++++- src/physics/cam/gw_rdg.F90 | 206 ++++++++++++++++++++++++++++++++++++ 2 files changed, 312 insertions(+), 4 deletions(-) diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index aeab27a5c6..627ec82e17 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -168,7 +168,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 :: & @@ -648,6 +650,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), & @@ -659,6 +663,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') @@ -734,15 +746,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', & @@ -1580,6 +1616,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. @@ -2257,6 +2299,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 @@ -2276,6 +2320,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) @@ -2306,6 +2351,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) @@ -2347,11 +2393,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 @@ -2385,6 +2432,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 @@ -2604,6 +2654,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) diff --git a/src/physics/cam/gw_rdg.F90 b/src/physics/cam/gw_rdg.F90 index b5a2a2137f..4c3e1ef745 100644 --- a/src/physics/cam/gw_rdg.F90 +++ b/src/physics/cam/gw_rdg.F90 @@ -19,6 +19,7 @@ module gw_rdg ! Public interface public :: gw_rdg_readnl public :: gw_rdg_src +public :: gw_rdg_resid_src public :: gw_rdg_belowpeak public :: gw_rdg_break_trap public :: gw_rdg_do_vdiff @@ -175,6 +176,211 @@ subroutine gw_rdg_readnl(nlfile) end subroutine gw_rdg_readnl +!========================================================================== +subroutine gw_rdg_resid_src(ncol, band, p, & + u, v, t, mxdis, kwvrdg, zi, nm, & + src_level, tend_level, tau, ubm, ubi, xv, yv, & + ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, c, tauoro ) + use gw_common, only: rair, GWBand + use gw_utils, only: dot_2d, midpoint_interp, get_unit_vector + !----------------------------------------------------------------------- + ! Orographic source for multiple gravity wave drag parameterization. + ! + ! The stress is returned for a single wave with c=0, over orography. + ! For points where the orographic variance is small (including ocean), + ! the returned stress is zero. + !------------------------------Arguments-------------------------------- + ! Column dimension. + integer, intent(in) :: ncol + + ! Band to emit orographic waves in. + ! Regardless, we will only ever emit into l = 0. + type(GWBand), intent(in) :: band + ! Pressure coordinates. + type(Coords1D), intent(in) :: p + + + ! Midpoint zonal/meridional winds. ( m s-1) + real(r8), intent(in) :: u(ncol,pver), v(ncol,pver) + ! Midpoint temperatures. (K) + real(r8), intent(in) :: t(ncol,pver) + ! Height estimate for ridge (m) [anisotropic orography]. + real(r8), intent(in) :: mxdis(ncol) + ! horiz wavenumber [anisotropic orography]. + real(r8), intent(in) :: kwvrdg(ncol) + ! Interface altitudes above ground (m). + real(r8), intent(in) :: zi(ncol,pver+1) + ! Midpoint Brunt-Vaisalla frequencies (s-1). + real(r8), intent(in) :: nm(ncol,pver) + + ! Indices of top gravity wave source level and lowest level where wind + ! tendencies are allowed. + integer, intent(out) :: src_level(ncol) + integer, intent(out) :: tend_level(ncol) + + ! Averages over source region. + real(r8), intent(out) :: nsrc(ncol) ! B-V frequency. + real(r8), intent(out) :: rsrc(ncol) ! Density. + real(r8), intent(out) :: usrc(ncol) ! Zonal wind. + real(r8), intent(out) :: vsrc(ncol) ! Meridional wind. + real(r8), intent(out) :: ubmsrc(ncol) ! On-obstacle wind. + ! normalized wavenumber + real(r8), intent(out) :: m2src(ncol) + + + ! Wave Reynolds stress. + real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) + ! Projection of wind at midpoints and interfaces. + real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1) + ! Unit vectors of source wind (zonal and meridional components). + real(r8), intent(out) :: xv(ncol), yv(ncol) + ! Phase speeds. + real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv) + ! source level mom. flux + real(r8), intent(out) :: tauoro(ncol) + + !---------------------------Local Storage------------------------------- + ! Column and level indices. + integer :: i, k + + ! Surface streamline displacement height (2*sgh). + real(r8) :: hdsp(ncol) + + ! Difference in interface pressure across source region. + real(r8) :: dpsrc(ncol) + ! Thickness of downslope wind region. + real(r8) :: ddw(ncol) + ! Thickness of linear wave region. + real(r8) :: dwv(ncol) + ! Wind speed in source region. + real(r8) :: wmsrc(ncol) + ! source level mom. flux + !real(r8) :: tauoro(ncol) + + real(r8) :: ragl(ncol) + real(r8) :: Fcrit_res,sghmax + +!-------------------------------------------------------------------------- +! Check that ngwav is equal to zero, otherwise end the job +!-------------------------------------------------------------------------- + if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not') + +!-------------------------------------------------------------------------- +! Average the basic state variables for the wave source over the depth of +! the orographic standard deviation. Here we assume that the appropiate +! values of wind, stability, etc. for determining the wave source are +! averages over the depth of the atmosphere penterated by the typical +! mountain. +! Reduces to the bottom midpoint values when mxdis=0, such as over ocean. +!-------------------------------------------------------------------------- + + Fcrit_res = 1.0_r8 + hdsp = mxdis ! no longer multipied by 2 + where(hdsp < 10._r8) + hdsp = 0._r8 + end where + + src_level = pver+1 + + tau(:,0,:) = 0.0_r8 + + ! Find depth of "source layer" for mountain waves + ! i.e., between ground and mountain top + do k = pver, 1, -1 + do i = 1, ncol + ! Need to have h >= z(k+1) here or code will bomb when h=0. + if ( (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k)) ) then + src_level(i) = k + end if + end do + end do + + rsrc = 0._r8 + usrc = 0._r8 + vsrc = 0._r8 + nsrc = 0._r8 + do i = 1, ncol + do k = pver, src_level(i), -1 + rsrc(i) = rsrc(i) + p%mid(i,k) / (rair*t(i,k))* p%del(i,k) + usrc(i) = usrc(i) + u(i,k) * p%del(i,k) + vsrc(i) = vsrc(i) + v(i,k) * p%del(i,k) + nsrc(i) = nsrc(i) + nm(i,k)* p%del(i,k) + end do + end do + + + do i = 1, ncol + dpsrc(i) = p%ifc(i,pver+1) - p%ifc(i,src_level(i)) + end do + + rsrc = rsrc / dpsrc + usrc = usrc / dpsrc + vsrc = vsrc / dpsrc + nsrc = nsrc / dpsrc + + ! Get the unit vector components and magnitude at the surface. + call get_unit_vector(usrc, vsrc, xv, yv, wmsrc ) + + ubmsrc = wmsrc + + ! Project the local wind at midpoints onto the source wind. + do k = 1, pver + ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv) + end do + + ! Compute the interface wind projection by averaging the midpoint winds. + ! Use the top level wind at the top interface. + ubi(:,1) = ubm(:,1) + + ubi(:,2:pver) = midpoint_interp(ubm) + + ! The minimum stratification allowing GW behavior + ! should really depend on horizontal scale since + ! + ! m^2 ~ (N/U)^2 - k^2 + ! + + m2src = ( (nsrc/(ubmsrc+0.01_r8))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01_r8))**2) + + ! Compute the interface wind projection by averaging the midpoint winds. + ! Use the top level wind at the top interface. + ubi(:,1) = ubm(:,1) + ubi(:,2:pver) = midpoint_interp(ubm) + ubi(:,pver+1) = ubm(:,pver) + + + + ! Determine the orographic c=0 source term following McFarlane (1987). + ! Set the source top interface index to pver, if the orographic term is + ! zero. + do i = 1, ncol + if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then + sghmax = Fcrit_res * (ubmsrc(i) / nsrc(i))**2 + tauoro(i) = 0.5_r8 * kwvrdg(i) * min(hdsp(i)**2, sghmax) * & + rsrc(i) * nsrc(i) * ubmsrc(i) + else + tauoro(i) = 0._r8 + end if + end do + + do i = 1, ncol + do k=src_level(i),pver+1 + tau(i,0,k) = tauoro(i) + end do + end do + + + ! Allow wind tendencies all the way to the model bottom. + tend_level = pver + + ! No spectrum; phase speed is just 0. + c = 0._r8 + +end subroutine gw_rdg_resid_src + + +!========================================================================== + subroutine gw_rdg_src(ncol, band, p, & u, v, t, mxdis, angxy, anixy, kwvrdg, iso, zi, nm, & src_level, tend_level, bwv_level ,tlb_level , tau, ubm, ubi, xv, yv, & From 176bbd3e22d243f32bfb5c2b7cf23a1e07c62521 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 23 Aug 2024 10:46:21 -0600 Subject: [PATCH 2/6] add template for two namelist variables for gw_drag --- bld/build-namelist | 5 +++++ bld/namelist_files/namelist_defaults_cam.xml | 4 ++++ bld/namelist_files/namelist_definition.xml | 12 +++++++++++ src/physics/cam/gw_drag.F90 | 22 +++++++++++++++----- 4 files changed, 38 insertions(+), 5 deletions(-) diff --git a/bld/build-namelist b/bld/build-namelist index bd91edc9f1..78e3c8c23c 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -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); @@ -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'); diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 308e128304..3374ad09b4 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -836,6 +836,10 @@ 0.5D0 0.5D0 +JULIO +JULIO +JULIO + 0.03D0 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index e0c4b46778..239c19f986 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -1338,6 +1338,12 @@ Whether or not to enable gravity waves from PBL moving mountains source. Default: .false. + +JULIO - need to add comment +Default: .false. + + Gravity wave spectrum dimension (wave numbers are from -pgwv to pgwv). @@ -1441,6 +1447,12 @@ Max efficiency associated with anisotropic OGW. Default: 1.0 + +JULIO - need to add comment and appropriate default +Default: ?????? JULIO + + Drag coefficient for obstacles in low-level flow. diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index 627ec82e17..3d50877cbe 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -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 @@ -252,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 @@ -361,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."// & @@ -746,7 +758,7 @@ 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' , & @@ -778,7 +790,7 @@ subroutine gw_init() 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', & @@ -1621,7 +1633,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) ! area fraction of res variance real(r8), pointer :: isowgt(:) - + ! Gamma ridges ! width of ridges. @@ -2659,7 +2671,7 @@ subroutine gw_rdg_calc( & 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, & From 9e2267cfbbf727d64d6f46118ab19bf5bae07430 Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Fri, 20 Dec 2024 14:36:08 -0700 Subject: [PATCH 3/6] commiting changes with residual orogrpahy. may not be 'perfect' --- bld/namelist_files/namelist_definition.xml | 2 +- src/physics/cam/gw_drag.F90 | 18 +++++++++++++----- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 239c19f986..6f07376438 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -1339,7 +1339,7 @@ Default: .false. + group="gw_drag_nl" valid_values="" > JULIO - need to add comment Default: .false. diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index 3d50877cbe..a71fac7141 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -111,9 +111,10 @@ module gw_drag 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 + ! Parameters controlling isotropic residual + ! orographic GW. + logical :: use_gw_rdg_resid = .false. + real(r8) :: effgw_rdg_resid = unset_r8 ! Horzontal wavelengths [m]. real(r8), parameter :: wavelength_mid = 1.e5_r8 @@ -2331,6 +2332,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) u, v, t, p, piln, zm, zi, & nm, ni, rhoi, kvtt, q, dse, & effgw_rdg_beta, effgw_rdg_beta_max, & + effgw_rdg_resid, use_gw_rdg_resid, & hwdth, clngt, gbxar, mxdis, angll, anixy, & isovar, isowgt, & rdg_beta_cd_llb, trpd_leewv_rdg_beta, & @@ -2362,6 +2364,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) u, v, t, p, piln, zm, zi, & nm, ni, rhoi, kvtt, q, dse, & effgw_rdg_gamma, effgw_rdg_gamma_max, & + effgw_rdg_resid, use_gw_rdg_resid, & hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, & isovar, isowgt, & rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, & @@ -2403,6 +2406,7 @@ subroutine gw_rdg_calc( & u, v, t, p, piln, zm, zi, & nm, ni, rhoi, kvtt, q, dse, & effgw_rdg, effgw_rdg_max, & + effgw_rdg_resid, luse_gw_rdg_resid, & hwdth, clngt, gbxar, & mxdis, angll, anixy, & isovar, isowgt, & @@ -2436,6 +2440,8 @@ subroutine gw_rdg_calc( & real(r8), intent(in) :: effgw_rdg ! Tendency efficiency. real(r8), intent(in) :: effgw_rdg_max + real(r8), intent(in) :: effgw_rdg_resid ! Tendency efficiency. + logical, intent(in) :: luse_gw_rdg_resid ! On-Off switch real(r8), intent(in) :: hwdth(ncol,prdg) ! width of ridges. real(r8), intent(in) :: clngt(ncol,prdg) ! length of ridges. real(r8), intent(in) :: gbxar(ncol) ! gridbox area @@ -2666,10 +2672,11 @@ subroutine gw_rdg_calc( & end do ! end of loop over multiple ridges + if (luse_gw_rdg_resid == .true.) then ! 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 + effgw = effgw_rdg_resid * isowgt !1.0_r8 * isowgt tauoro = 0._r8 call gw_rdg_resid_src(ncol, band_oro, p, & @@ -2716,7 +2723,8 @@ subroutine gw_rdg_calc( & 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 + ! end of residual variance calc + end if ! Calculate energy change for output to CAM's energy checker. call energy_change(dt, p, u, v, ptend%u(:ncol,:), & From 3068bad2f9718eee6563d46212eb8b86e8d6456d Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Mon, 23 Dec 2024 11:30:34 -0700 Subject: [PATCH 4/6] minor mods to gw_drag --- src/physics/cam/gw_drag.F90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index a71fac7141..74274ea0af 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -812,6 +812,12 @@ subroutine gw_init() 'Ridge based momentum flux profile') call register_vector_field('TAUARDGBETAX','TAUARDGBETAY') + call addfld('TAURESIDBETAY' , (/ 'ilev' /) , 'I' ,'N m-2' , & + 'Ridge based momentum flux profile') + call addfld('TAURESIDBETAX' , (/ 'ilev' /) , 'I' ,'N m-2' , & + 'Ridge based momentum flux profile') + call register_vector_field('TAURESIDBETAX','TAURESIDBETAY') + if (history_waccm) then call add_default('TAUARDGBETAX', 1, ' ') call add_default('TAUARDGBETAY ', 1, ' ') @@ -2672,6 +2678,9 @@ subroutine gw_rdg_calc( & end do ! end of loop over multiple ridges + call outfld('TAUARDG'//trim(type)//'X', taurx, ncol, lchnk) + call outfld('TAUARDG'//trim(type)//'Y', taury, ncol, lchnk) + if (luse_gw_rdg_resid == .true.) then ! Add additional GW from residual variance. Assumed isotropic !kwvrdg = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!! @@ -2720,6 +2729,9 @@ subroutine gw_rdg_calc( & 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('TAURESID'//trim(type)//'X', taurx, ncol, lchnk) + call outfld('TAURESID'//trim(type)//'Y', taury, 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) @@ -2731,8 +2743,6 @@ subroutine gw_rdg_calc( & ptend%v(:ncol,:), ptend%s(:ncol,:), de) flx_heat(:ncol) = de - call outfld('TAUARDG'//trim(type)//'X', taurx, ncol, lchnk) - call outfld('TAUARDG'//trim(type)//'Y', taury, ncol, lchnk) if (trim(type) == 'BETA') then fname(1) = 'TAUGWX' From 04f41af4f525fe17551fa3394ac1aad969e4f7bc Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Thu, 9 Jan 2025 11:44:48 -0700 Subject: [PATCH 5/6] added vorticity coupling from SE dycore to GW param --- src/dynamics/se/dp_coupling.F90 | 33 ++++++- src/dynamics/se/dyn_comp.F90 | 4 + src/dynamics/se/gravity_waves_sources.F90 | 107 +++++++++++++++++++++- src/physics/cam/gw_common.F90 | 5 +- src/physics/cam/gw_drag.F90 | 27 +++++- src/physics/cam/gw_movmtn.F90 | 69 +++++++++++--- 6 files changed, 221 insertions(+), 24 deletions(-) diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 919b7f3510..e7b9331083 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -50,7 +50,7 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) ! dry air mass. use gravity_waves_sources, only: gws_src_fnct - use dyn_comp, only: frontgf_idx, frontga_idx + use dyn_comp, only: frontgf_idx, frontga_idx, vort4gw_idx use phys_control, only: use_gw_front, use_gw_front_igw use hycoef, only: hyai, ps0 use fvm_mapping, only: dyn2phys_vector, dyn2phys_all_vars @@ -84,9 +84,18 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) real (kind=r8), allocatable :: frontga(:,:,:) ! function (frontgf) and angle (frontga) real (kind=r8), allocatable :: frontgf_phys(:,:,:) real (kind=r8), allocatable :: frontga_phys(:,:,:) + + !++jtb 12/31/24 + ! Vorticity + real (kind=r8), allocatable :: vort4gw(:,:,:) ! temp arrays to hold vorticity + real (kind=r8), allocatable :: vort4gw_phys(:,:,:) + + ! Pointers to pbuf real (kind=r8), pointer :: pbuf_frontgf(:,:) real (kind=r8), pointer :: pbuf_frontga(:,:) + !++jtb 12/31/24 + real (kind=r8), pointer :: pbuf_vort4gw(:,:) integer :: ncols, ierr integer :: col_ind, blk_ind(1), m @@ -110,6 +119,10 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) nullify(pbuf_chnk) nullify(pbuf_frontgf) nullify(pbuf_frontga) + !++jtb + nullify(pbuf_vort4gw) + + if (fv_nphys > 0) then nphys = fv_nphys @@ -135,11 +148,14 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (ierr /= 0) call endrun("dp_coupling: Allocate of frontgf failed.") allocate(frontga(nphys_pts,pver,nelemd), stat=ierr) if (ierr /= 0) call endrun("dp_coupling: Allocate of frontga failed.") + !++jtb + allocate(vort4gw(nphys_pts,pver,nelemd), stat=ierr) + if (ierr /= 0) call endrun("dp_coupling: Allocate of vort4gw failed.") end if if (iam < par%nprocs) then if (use_gw_front .or. use_gw_front_igw) then - call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys) + call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, vort4gw, nphys) end if if (fv_nphys > 0) then @@ -204,6 +220,8 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then frontgf(:,:,:) = 0._r8 frontga(:,:,:) = 0._r8 + !++jtb + vort4gw(:,:,:) = 0._r8 end if endif ! iam < par%nprocs @@ -222,6 +240,8 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then allocate(frontgf_phys(pcols, pver, begchunk:endchunk)) allocate(frontga_phys(pcols, pver, begchunk:endchunk)) + !++jtb 12/31/24 + allocate(vort4gw_phys(pcols, pver, begchunk:endchunk)) end if !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m) do col_ind = 1, phys_columns_on_task @@ -239,6 +259,8 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then frontgf_phys(icol, ilyr, lchnk) = frontgf(blk_ind(1), ilyr, ie) frontga_phys(icol, ilyr, lchnk) = frontga(blk_ind(1), ilyr, ie) + !++jtb 12/31/24 + vort4gw_phys(icol, ilyr, lchnk) = vort4gw(blk_ind(1), ilyr, ie) end if end do @@ -248,22 +270,27 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) end do end do end do + !++jtb 12/31/24 if (use_gw_front .or. use_gw_front_igw) then - !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga) + !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga, pbuf_vort4gw) do lchnk = begchunk, endchunk ncols = get_ncols_p(lchnk) pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf) call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga) + call pbuf_get_field(pbuf_chnk, vort4gw_idx, pbuf_vort4gw) do icol = 1, ncols do ilyr = 1, pver pbuf_frontgf(icol, ilyr) = frontgf_phys(icol, ilyr, lchnk) pbuf_frontga(icol, ilyr) = frontga_phys(icol, ilyr, lchnk) + pbuf_vort4gw(icol, ilyr) = vort4gw_phys(icol, ilyr, lchnk) end do end do end do deallocate(frontgf_phys) deallocate(frontga_phys) + !++jtb 12/31/24 + deallocate(vort4gw_phys) end if call t_stopf('dpcopy') diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 586ee06b1f..56f7ee5355 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -79,6 +79,8 @@ module dyn_comp ! Frontogenesis indices integer, public :: frontgf_idx = -1 integer, public :: frontga_idx = -1 +!++jtb +integer, public :: vort4gw_idx = -1 interface read_dyn_var module procedure read_dyn_field_2d @@ -571,6 +573,8 @@ subroutine dyn_register() frontgf_idx) call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/), & frontga_idx) + call pbuf_add_field("VORT4GW", "global", dtype_r8, (/pcols,pver/), & + vort4gw_idx) end if end subroutine dyn_register diff --git a/src/dynamics/se/gravity_waves_sources.F90 b/src/dynamics/se/gravity_waves_sources.F90 index a929dfeaf1..2789a7e385 100644 --- a/src/dynamics/se/gravity_waves_sources.F90 +++ b/src/dynamics/se/gravity_waves_sources.F90 @@ -6,6 +6,9 @@ module gravity_waves_sources use hybrid_mod, only: hybrid_t use shr_kind_mod, only: r8 => shr_kind_r8 + !++ jtb (added for now, while debugging) + use cam_logfile, only: iulog + implicit none private save @@ -18,6 +21,7 @@ module gravity_waves_sources public :: gws_src_fnct public :: gws_init private :: compute_frontogenesis + private :: compute_vorticity_4gw type (EdgeBuffer_t) :: edge3 type (derivative_t) :: deriv @@ -45,7 +49,7 @@ subroutine gws_init(elem) end subroutine gws_init - subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) + subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga, vort4gw, nphys) use derivative_mod, only : derivinit use dimensions_mod, only : npsq, nelemd use dof_mod, only : UniquePoints @@ -60,11 +64,21 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) real (kind=r8), intent(out) :: frontgf(nphys*nphys,pver,nelemd) real (kind=r8), intent(out) :: frontga(nphys*nphys,pver,nelemd) + !++jtb (12/31/24) + real (kind=r8), intent(out) :: vort4gw(nphys*nphys,pver,nelemd) + !!real (kind=r8) :: vort4gw(nphys*nphys,pver,nelemd) + + ! Local variables type (hybrid_t) :: hybrid integer :: nets, nete, ithr, ncols, ie real(kind=r8), allocatable :: frontgf_thr(:,:,:,:) real(kind=r8), allocatable :: frontga_thr(:,:,:,:) + + !++jtb (12/31/24) + real(kind=r8), allocatable :: vort4gw_thr(:,:,:,:) + + ! This does not need to be a thread private data-structure call derivinit(deriv) @@ -75,25 +89,114 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) allocate(frontgf_thr(nphys,nphys,nlev,nets:nete)) allocate(frontga_thr(nphys,nphys,nlev,nets:nete)) + !++jtb (12/31/24) + allocate(vort4gw_thr(nphys,nphys,nlev,nets:nete)) + call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys) + call compute_vorticity_4gw(vort4gw_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys) + if (fv_nphys>0) then do ie=nets,nete frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/)) frontga(:,:,ie) = RESHAPE(frontga_thr(:,:,:,ie),(/nphys*nphys,nlev/)) + !++jtb (12/31/24) + vort4gw(:,:,ie) = RESHAPE(vort4gw_thr(:,:,:,ie),(/nphys*nphys,nlev/)) end do else do ie=nets,nete ncols = elem(ie)%idxP%NumUniquePts call UniquePoints(elem(ie)%idxP, nlev, frontgf_thr(:,:,:,ie), frontgf(1:ncols,:,ie)) call UniquePoints(elem(ie)%idxP, nlev, frontga_thr(:,:,:,ie), frontga(1:ncols,:,ie)) + !++jtb (12/31/24) + call UniquePoints(elem(ie)%idxP, nlev, vort4gw_thr(:,:,:,ie), vort4gw(1:ncols,:,ie)) end do end if deallocate(frontga_thr) deallocate(frontgf_thr) + !++ jtb 12/31/24 + deallocate(vort4gw_thr) + !!$OMP END PARALLEL end subroutine gws_src_fnct + !++jtb (12/31/24) + subroutine compute_vorticity_4gw(vort4gw,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! compute vorticity for use in gw params + ! F = ( curl ) [U,V] + ! + ! Original by Peter Lauritzen, Julio Bacmeister*, Dec 2024 + ! Patterned on 'compute_frontogenesis' + ! + ! * corresponding/blame-able + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + use derivative_mod, only: vorticity_sphere + use edge_mod, only: edgevpack, edgevunpack + use bndry_mod, only: bndry_exchange + use dyn_grid, only: hvcoord + use dimensions_mod, only: fv_nphys,ntrac + use fvm_mapping, only: dyn2phys_vector,dyn2phys + + type(hybrid_t), intent(in) :: hybrid + type(element_t), intent(inout), target :: elem(:) + type(derivative_t), intent(in) :: ederiv + integer, intent(in) :: nets,nete,nphys + integer, intent(in) :: tl,tlq + real(r8), intent(out) :: vort4gw(nphys,nphys,nlev,nets:nete) + + ! local + real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np) + !!real(r8) :: vort_tmp(fv_nphys*fv_nphys,nlev) + real(r8) :: vort_gll(np,np,nlev,nets:nete) + integer :: k,kptr,i,j,ie,component,h,nq,m_cnst,n0 + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! First calculate vorticity on GLL grid + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! set timelevel=1 fro velocities + n0=1 + do ie=nets,nete + do k=1,nlev + call vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),ederiv,elem(ie),vort_gll(:,:,k,ie)) + end do + do k=1,nlev + vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%spheremp(:,:) + end do + ! pack ++jtb no idea what these routines are doing + call edgeVpack(edge3, vort_gll(:,:,:,ie),nlev,0,ie) + enddo + call bndry_exchange(hybrid,edge3,location='compute_vorticity_4gw') + do ie=nets,nete + call edgeVunpack(edge3, vort_gll(:,:,:,ie),nlev,0,ie) + ! apply inverse mass matrix, + do k=1,nlev + vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:) + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Now regrid from GLL to PhysGrid if necessary + ! otherwise just return vorticity on GLL grid + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (fv_nphys>0) then + tmp = 1.0_r8 + area_inv = dyn2phys(tmp,elem(ie)%metdet) + area_inv = 1.0_r8/area_inv + !!! vort_tmp(:,:) = dyn2phys(vort_gll(:,:,:,ie),elem(ie)) !peter replace with scalar mapping !++jtb: Think I did that ... + do k=1,nlev + vort4gw(:,:,k,ie) = dyn2phys( vort_gll(:,:,k,ie) , elem(ie)%metdet , area_inv ) + end do + else + do k=1,nlev + vort4gw(:,:,k,ie) = vort_gll(:,:,k,ie) + end do + end if + enddo + + + end subroutine compute_vorticity_4gw + + subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! compute frontogenesis function F @@ -170,6 +273,8 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets, enddo ! pack call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie) + !++jtb: + ! Why are dims 2*nlev,nlev, not 2*nlev,2*nlev, or nlev,nlev, ???? call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie) enddo call bndry_exchange(hybrid,edge3,location='compute_frontogenesis') diff --git a/src/physics/cam/gw_common.F90 b/src/physics/cam/gw_common.F90 index 04014c8c97..b98d8155b8 100644 --- a/src/physics/cam/gw_common.F90 +++ b/src/physics/cam/gw_common.F90 @@ -132,7 +132,10 @@ function new_GWBand(ngwv, dc, fcrit2, wavelength) result(band) ! Simple assignments. band%ngwv = ngwv band%dc = dc - band%fcrit2 = fcrit2 + !++jtb (12/24/2024) + ! This nonsesnse needs to be straightened out + ! For now just ensure fcrit is always set to 1 + band%fcrit2 = 1.0_r8 ! fcrit2 ! Uniform phase speed reference grid. allocate(band%cref(-ngwv:ngwv)) diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index 212051bcde..d5192c8e77 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -161,6 +161,9 @@ module gw_drag integer :: ttend_sh_idx = -1 integer :: frontgf_idx = -1 integer :: frontga_idx = -1 + !++jtb + integer :: vort4gw_idx = -1 + integer :: sgh_idx = -1 ! From CLUBB @@ -367,10 +370,12 @@ subroutine gw_drag_readnl(nlfile) call mpi_bcast(effgw_rdg_resid, 1, mpi_real8, mstrid, mpicom, ierr) if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_resid") + !++ jtb (12/24/2024) + ! This is confusing. Not sure when if ever this was needed. ! Check if fcrit2 was set. - call shr_assert(fcrit2 /= unset_r8, & - "gw_drag_readnl: fcrit2 must be set via the namelist."// & - errMsg(__FILE__, __LINE__)) + !call shr_assert(fcrit2 /= unset_r8, & + ! "gw_drag_readnl: fcrit2 must be set via the namelist."// & + ! errMsg(__FILE__, __LINE__)) ! Check if pgwv was set. call shr_assert(pgwv >= 0, & @@ -911,6 +916,7 @@ subroutine gw_init() frontgf_idx = pbuf_get_index('FRONTGF') frontga_idx = pbuf_get_index('FRONTGA') + vort4gw_idx = pbuf_get_index('VORT4GW') call shr_assert(unset_r8 /= frontgfc, & "gw_init: Frontogenesis enabled, but frontgfc was & @@ -935,6 +941,8 @@ subroutine gw_init() 'Frontogenesis function at gws src level') call addfld ('FRONTGFA', (/ 'lev' /), 'A', 'K^2/M^2/S', & 'Frontogenesis function at gws src level') + call addfld ('VORT4GW', (/ 'lev' /), 'A', '1/S', & + 'Vorticity') if (history_waccm) then call add_default('FRONTGF', 1, ' ') @@ -1596,6 +1604,8 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) ! Frontogenesis real(r8), pointer :: frontgf(:,:) real(r8), pointer :: frontga(:,:) + !++jtb 12/31/24 + real(r8), pointer :: vort4gw(:,:) ! Temperature change due to deep convection. real(r8), pointer :: ttend_dp(:,:) @@ -1785,11 +1795,17 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) call pbuf_get_field(pbuf, upwp_clubb_gw_idx, upwp_clubb_gw) call pbuf_get_field(pbuf, vpwp_clubb_gw_idx, vpwp_clubb_gw) + !++jtb 01/03/25 + ! Vorticity from SE dycore. This needs to be either + ! generalized to other dycores or protected with some + ! endrun if dycore != SE + call pbuf_get_field(pbuf, vort4gw_idx, vort4gw) + xpwp_clubb(:ncol,:) = sqrt( upwp_clubb_gw(:ncol,:)**2 + vpwp_clubb_gw(:ncol,:)**2 ) effgw = 1._r8 call gw_movmtn_src(ncol, lchnk, band_movmtn , movmtn_desc, & - u, v, ttend_dp(:ncol,:), ttend_clubb(:ncol,:), xpwp_clubb(:ncol,:) , & + u, v, ttend_dp(:ncol,:), ttend_clubb(:ncol,:), xpwp_clubb(:ncol,:), vort4gw(:ncol,:), & zm, alpha_gw_movmtn, src_level, tend_level, & tau, ubm, ubi, xv, yv, & phase_speeds, hdepth) @@ -1849,6 +1865,9 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) call outfld('UPWP_CLUBB_GW', upwp_clubb_gw, pcols, lchnk) call outfld('VPWP_CLUBB_GW', vpwp_clubb_gw, pcols, lchnk) + !++jtb 01/03/25 (see comment above) + call outfld ('VORT4GW', vort4gw, pcols, lchnk) + !Deallocate variables that are no longer used: deallocate(tau, gwut, phase_speeds) diff --git a/src/physics/cam/gw_movmtn.F90 b/src/physics/cam/gw_movmtn.F90 index 0408928932..27d04983a2 100644 --- a/src/physics/cam/gw_movmtn.F90 +++ b/src/physics/cam/gw_movmtn.F90 @@ -2,7 +2,7 @@ module gw_movmtn ! ! This module parameterizes gravity waves generated by the obstacle effect produced by -! boundary layer turbulence for convection. +! internal circulations in the atmosphere. ! use gw_utils, only: r8 @@ -35,12 +35,12 @@ module gw_movmtn !========================================================================== subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, & - netdt, netdt_shcu, xpwp_shcu, & + netdt, netdt_shcu, xpwp_shcu, vorticity, & zm, alpha_gw_movmtn, src_level, tend_level, tau, ubm, ubi, xv, yv, & c, hdepth) !----------------------------------------------------------------------- ! Flexible driver for gravity wave source from obstacle effects produced -! by boundary layer turbulence or deep convection +! by internal circulations !----------------------------------------------------------------------- use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp use gw_common, only: GWBand, pver, qbo_hdepth_scaling @@ -65,6 +65,8 @@ subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, & real(r8), intent(in) :: netdt_shcu(:,:) ! Higher order flux from ShCu/PBL. real(r8), intent(in) :: xpwp_shcu(ncol,pver+1) + ! Relative vorticity + real(r8), intent(in) :: vorticity(ncol,pver) ! Midpoint altitudes. real(r8), intent(in) :: zm(ncol,pver) ! tunable parameter controlling proportion of PBL momentum flux emitted as GW @@ -136,10 +138,12 @@ subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, & ! Index for ground based phase speed bin real(r8) :: c0(ncol,-band%ngwv:band%ngwv) integer :: c_idx(ncol,-band%ngwv:band%ngwv) - ! Flux source from ShCu/PBL + ! GW Flux source from ShCu/PBL real(r8) :: xpwp_src(ncol) + ! Gw Flux source from vorticity + real(r8) :: vort_src(ncol) ! Manual steering level set - integer :: Steer_k + integer :: Steer_k, Steer_k_vort !---------------------------------------------------------------------- ! Initialize tau array @@ -150,18 +154,21 @@ subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, & tau0 = 0.0_r8 !---------------------------------------------------------------------- - ! Calculate flux source from ShCu/PBL + ! Calculate flux source from ShCu/PBL and set Steering level !---------------------------------------------------------------------- - xpwp_src = shcu_flux_src( xpwp_shcu, ncol, pver+1, alpha_gw_movmtn ) - + ! call shcu_flux_src( xpwp_shcu, ncol, pver+1, alpha_gw_movmtn, xpwp_src, Steer_k ) + + !---------------------------------------------------------------------- + ! Calculate flux source from vorticity + !---------------------------------------------------------------------- + call vorticity_flux_src( vorticity, ncol, pver , alpha_gw_movmtn, xpwp_src, Steer_k ) + !------------------------------------------------------------------------ ! Determine wind and unit vectors approximately at the source (steering level), then ! project winds. !------------------------------------------------------------------------ - ! Winds at 'steering level' - Steer_k = pver-1 - usteer = u(:,Steer_k) !k defined in line21 (at specified altitude) + usteer = u(:,Steer_k) ! vsteer = v(:,Steer_k) steer_level = real(Steer_k,r8) @@ -209,7 +216,7 @@ subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, & if (use_gw_movmtn_pbl) then boti=pver - topi=Steer_k-10 ! desc%k-5 + topi=Steer_k-10 !++ tuning test 01/05/25 else do k = pver, 1, -1 !start at surface do i = 1, ncol @@ -419,15 +426,19 @@ pure function index_of_nearest(x, grid) result(idx) end function index_of_nearest !!!!!!!!!!!!!!!!!!!!!!!!!!! -pure function shcu_flux_src (xpwp_shcu , ncol, pverx, alpha_gw_movmtn ) result(xpwp_src) +subroutine shcu_flux_src (xpwp_shcu , ncol, pverx, alpha_gw_movmtn, xpwp_src, steering_level ) !! result(xpwp_src) + !!! use gw_common, only: pver integer, intent(in) :: ncol,pverx real(r8), intent(in) :: xpwp_shcu (ncol,pverx) real(r8), intent(in) :: alpha_gw_movmtn - real(r8) :: xpwp_src(ncol) + real(r8), intent(out) :: xpwp_src(ncol) + integer, intent(out) :: steering_level integer :: k, nlayers + steering_level = (pverx-1) - 5 !++ tuning test 12/30/24 + !----------------------------------- ! Simple average over layers. ! Probably can do better @@ -439,6 +450,34 @@ pure function shcu_flux_src (xpwp_shcu , ncol, pverx, alpha_gw_movmtn ) result(x end do xpwp_src(:) = alpha_gw_movmtn * xpwp_src(:)/(1.0_r8*nlayers) -end function shcu_flux_src +end subroutine shcu_flux_src + +!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine vorticity_flux_src (vorticity , ncol, pverx, alpha_gw_movmtn, vort_src, steering_level ) !! result(xpwp_src) + !!! use gw_common, only: pver + integer, intent(in) :: ncol,pverx + real(r8), intent(in) :: vorticity (ncol,pverx) + real(r8), intent(in) :: alpha_gw_movmtn + + real(r8), intent(out) :: vort_src(ncol) + integer, intent(out) :: steering_level + + real(r8) :: scale_factor + integer :: k, nlayers + + steering_level = pverx - 10 !++ ????? + scale_factor = 1.e4 ! scales vorticity amp to u'w' in CLUBB + !----------------------------------- + ! Simple average over layers. + ! Probably can do better + !----------------------------------- + nlayers=10 + vort_src(:) =0._r8 + do k = 0, nlayers-1 + vort_src(:) = vort_src(:) + scale_factor * abs( vorticity(:,pverx-k) ) + end do + vort_src(:) = alpha_gw_movmtn * vort_src(:)/(1.0_r8*nlayers) + +end subroutine vorticity_flux_src end module gw_movmtn From f9986872ff98853dc1aab10053315566072fa36a Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Wed, 15 Jan 2025 07:35:01 -0700 Subject: [PATCH 6/6] Latest movmtn mods - vorticity source --- src/dynamics/se/dp_coupling.F90 | 42 +++++++++++++++++------ src/dynamics/se/dyn_comp.F90 | 8 +++-- src/dynamics/se/gravity_waves_sources.F90 | 11 +++--- src/physics/cam/gw_drag.F90 | 8 +++-- src/physics/cam/gw_movmtn.F90 | 2 +- 5 files changed, 49 insertions(+), 22 deletions(-) diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index e7b9331083..2b339ca80e 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -51,7 +51,7 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) use gravity_waves_sources, only: gws_src_fnct use dyn_comp, only: frontgf_idx, frontga_idx, vort4gw_idx - use phys_control, only: use_gw_front, use_gw_front_igw + use phys_control, only: use_gw_front, use_gw_front_igw, use_gw_movmtn_pbl use hycoef, only: hyai, ps0 use fvm_mapping, only: dyn2phys_vector, dyn2phys_all_vars use se_dyn_time_mod, only: timelevel_qdp @@ -85,7 +85,7 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) real (kind=r8), allocatable :: frontgf_phys(:,:,:) real (kind=r8), allocatable :: frontga_phys(:,:,:) - !++jtb 12/31/24 + !++jtb 01/14/25 ! Vorticity real (kind=r8), allocatable :: vort4gw(:,:,:) ! temp arrays to hold vorticity real (kind=r8), allocatable :: vort4gw_phys(:,:,:) @@ -148,13 +148,16 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (ierr /= 0) call endrun("dp_coupling: Allocate of frontgf failed.") allocate(frontga(nphys_pts,pver,nelemd), stat=ierr) if (ierr /= 0) call endrun("dp_coupling: Allocate of frontga failed.") - !++jtb + end if + if (use_gw_movmtn_pbl) then + !++jtb 01/14/25 allocate(vort4gw(nphys_pts,pver,nelemd), stat=ierr) if (ierr /= 0) call endrun("dp_coupling: Allocate of vort4gw failed.") end if + !++jtb 01/14/25 if (iam < par%nprocs) then - if (use_gw_front .or. use_gw_front_igw) then + if (use_gw_front .or. use_gw_front_igw .or. use_gw_movmtn_pbl ) then call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, vort4gw, nphys) end if @@ -220,7 +223,9 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then frontgf(:,:,:) = 0._r8 frontga(:,:,:) = 0._r8 - !++jtb + end if + if (use_gw_movmtn_pbl) then + !++jtb 01/14/25 vort4gw(:,:,:) = 0._r8 end if @@ -240,7 +245,9 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then allocate(frontgf_phys(pcols, pver, begchunk:endchunk)) allocate(frontga_phys(pcols, pver, begchunk:endchunk)) - !++jtb 12/31/24 + end if + if (use_gw_movmtn_pbl) then + !++jtb 01/14/25 allocate(vort4gw_phys(pcols, pver, begchunk:endchunk)) end if !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m) @@ -259,7 +266,9 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) if (use_gw_front .or. use_gw_front_igw) then frontgf_phys(icol, ilyr, lchnk) = frontgf(blk_ind(1), ilyr, ie) frontga_phys(icol, ilyr, lchnk) = frontga(blk_ind(1), ilyr, ie) - !++jtb 12/31/24 + end if + if (use_gw_movmtn_pbl) then + !++jtb 01/14/25 vort4gw_phys(icol, ilyr, lchnk) = vort4gw(blk_ind(1), ilyr, ie) end if end do @@ -270,7 +279,6 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) end do end do end do - !++jtb 12/31/24 if (use_gw_front .or. use_gw_front_igw) then !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga, pbuf_vort4gw) do lchnk = begchunk, endchunk @@ -278,18 +286,30 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf) call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga) - call pbuf_get_field(pbuf_chnk, vort4gw_idx, pbuf_vort4gw) do icol = 1, ncols do ilyr = 1, pver pbuf_frontgf(icol, ilyr) = frontgf_phys(icol, ilyr, lchnk) pbuf_frontga(icol, ilyr) = frontga_phys(icol, ilyr, lchnk) - pbuf_vort4gw(icol, ilyr) = vort4gw_phys(icol, ilyr, lchnk) end do end do end do deallocate(frontgf_phys) deallocate(frontga_phys) - !++jtb 12/31/24 + end if + !++jtb 01/14/25 + if (use_gw_movmtn_pbl) then + !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga, pbuf_vort4gw) + do lchnk = begchunk, endchunk + ncols = get_ncols_p(lchnk) + pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) + call pbuf_get_field(pbuf_chnk, vort4gw_idx, pbuf_vort4gw) + do icol = 1, ncols + do ilyr = 1, pver + pbuf_vort4gw(icol, ilyr) = vort4gw_phys(icol, ilyr, lchnk) + end do + end do + end do + !++jtb 01/14/25 deallocate(vort4gw_phys) end if diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 56f7ee5355..d5e22f8ce9 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -10,7 +10,7 @@ module dyn_comp cnst_is_a_water_species use cam_control_mod, only: initial_run use cam_initfiles, only: initial_file_get_id, topo_file_get_id, pertlim -use phys_control, only: use_gw_front, use_gw_front_igw +use phys_control, only: use_gw_front, use_gw_front_igw, use_gw_movmtn_pbl use dyn_grid, only: ini_grid_name, timelevel, hvcoord, edgebuf, & ini_grid_hdim_name @@ -573,6 +573,8 @@ subroutine dyn_register() frontgf_idx) call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/), & frontga_idx) + end if + if (use_gw_movmtn_pbl) then call pbuf_add_field("VORT4GW", "global", dtype_r8, (/pcols,pver/), & vort4gw_idx) end if @@ -879,8 +881,8 @@ subroutine dyn_init(dyn_in, dyn_out) call get_loop_ranges(hybrid, ibeg=nets, iend=nete) call prim_init2(elem, fvm, hybrid, nets, nete, TimeLevel, hvcoord) !$OMP END PARALLEL - - if (use_gw_front .or. use_gw_front_igw) call gws_init(elem) + !++jtb 01/14/25 + if (use_gw_front .or. use_gw_front_igw .or. use_gw_movmtn_pbl) call gws_init(elem) end if ! iam < par%nprocs call addfld ('nu_kmvis', (/ 'lev' /), 'A', '', 'Molecular viscosity Laplacian coefficient' , gridname='GLL') diff --git a/src/dynamics/se/gravity_waves_sources.F90 b/src/dynamics/se/gravity_waves_sources.F90 index 2789a7e385..2d90b759b2 100644 --- a/src/dynamics/se/gravity_waves_sources.F90 +++ b/src/dynamics/se/gravity_waves_sources.F90 @@ -23,7 +23,7 @@ module gravity_waves_sources private :: compute_frontogenesis private :: compute_vorticity_4gw - type (EdgeBuffer_t) :: edge3 + type (EdgeBuffer_t) :: edge3,edge1 type (derivative_t) :: deriv real(r8) :: psurf_ref @@ -44,6 +44,7 @@ subroutine gws_init(elem) ! Set up variables similar to dyn_comp and prim_driver_mod initializations call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1) + call initEdgeBuffer(par, edge1, elem, nlev,nthreads=1) psurf_ref = hypi(plev+1) @@ -66,7 +67,7 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga, vort4gw, nphys) !++jtb (12/31/24) real (kind=r8), intent(out) :: vort4gw(nphys*nphys,pver,nelemd) - !!real (kind=r8) :: vort4gw(nphys*nphys,pver,nelemd) + !!real (kind=r8) :: vort4gw(nphys*nphys,pver,nelemd) phl remove ! Local variables @@ -164,11 +165,11 @@ subroutine compute_vorticity_4gw(vort4gw,tl,tlq,elem,ederiv,hybrid,nets,nete,nph vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%spheremp(:,:) end do ! pack ++jtb no idea what these routines are doing - call edgeVpack(edge3, vort_gll(:,:,:,ie),nlev,0,ie) + call edgeVpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie) enddo - call bndry_exchange(hybrid,edge3,location='compute_vorticity_4gw') + call bndry_exchange(hybrid,edge1,location='compute_vorticity_4gw') do ie=nets,nete - call edgeVunpack(edge3, vort_gll(:,:,:,ie),nlev,0,ie) + call edgeVunpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie) ! apply inverse mass matrix, do k=1,nlev vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:) diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90 index d5192c8e77..996b47c5a8 100644 --- a/src/physics/cam/gw_drag.F90 +++ b/src/physics/cam/gw_drag.F90 @@ -1534,7 +1534,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) use gw_front, only: gw_cm_src use gw_convect, only: gw_beres_src use gw_movmtn, only: gw_movmtn_src - + use dycore, only: dycore_is !------------------------------Arguments-------------------------------- type(physics_state), intent(in) :: state ! physics state structure type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer @@ -1799,7 +1799,11 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat) ! Vorticity from SE dycore. This needs to be either ! generalized to other dycores or protected with some ! endrun if dycore != SE - call pbuf_get_field(pbuf, vort4gw_idx, vort4gw) + if (dycore_is('SE')) then + call pbuf_get_field(pbuf, vort4gw_idx, vort4gw) + else + call endrun( 'gw_drag: vort4gw only with SE') + end if xpwp_clubb(:ncol,:) = sqrt( upwp_clubb_gw(:ncol,:)**2 + vpwp_clubb_gw(:ncol,:)**2 ) diff --git a/src/physics/cam/gw_movmtn.F90 b/src/physics/cam/gw_movmtn.F90 index 27d04983a2..e4854b23c4 100644 --- a/src/physics/cam/gw_movmtn.F90 +++ b/src/physics/cam/gw_movmtn.F90 @@ -465,7 +465,7 @@ subroutine vorticity_flux_src (vorticity , ncol, pverx, alpha_gw_movmtn, vort_sr real(r8) :: scale_factor integer :: k, nlayers - steering_level = pverx - 10 !++ ????? + steering_level = pverx - 20 !++ ????? scale_factor = 1.e4 ! scales vorticity amp to u'w' in CLUBB !----------------------------------- ! Simple average over layers.