diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist index 8dc532b6a3df..670390c6f66e 100755 --- a/components/eam/bld/build-namelist +++ b/components/eam/bld/build-namelist @@ -830,6 +830,9 @@ add_default($nl,'use_hetfrz_classnuc'); add_default($nl,'hist_hetfrz_classnuc'); add_default($nl,'gw_convect_hcf') if (get_default_value('gw_convect_hcf')); add_default($nl,'hdepth_scaling_factor') if (get_default_value('hdepth_scaling_factor')); +add_default($nl,'gw_convect_hdepth_min') if (get_default_value('gw_convect_hdepth_min')); +add_default($nl,'gw_convect_storm_speed_min') if (get_default_value('gw_convect_storm_speed_min')); +add_default($nl,'use_gw_convect_old', 'val'=>'.true.'); add_default($nl,'linoz_psc_T'); if ($cfg->get('microphys') =~ /^mg2/) { add_default($nl,'micro_mg_dcs_tdep'); diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index 44e4e7e65ae5..ba4b30bf4465 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -1890,6 +1890,8 @@ with se_tstep, dt_remap_factor, dt_tracer_factor set to -1 10.0 0.50 1.0 +2.5 +10.0 0.375 .true. 2.5D0 diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index 8228c7d8d2e3..6c00afd21e50 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -1110,6 +1110,23 @@ Scaling factor for the heating depth Default: 1.0 + +minimum hdepth for for convective GWD spectrum lookup table +Default: 2.5 + + + +minimum convective storm speed for convective GWD +Default: 10.0 m/s + + + +switch to revert to old calculation of Beres scheme for heating depth and max + + Efficiency associated with convective gravity waves from the Beres diff --git a/components/eam/src/physics/cam/gw_convect.F90 b/components/eam/src/physics/cam/gw_convect.F90 index fc2fce99214a..82d6417ce811 100644 --- a/components/eam/src/physics/cam/gw_convect.F90 +++ b/components/eam/src/physics/cam/gw_convect.F90 @@ -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 @@ -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(:,:,:) @@ -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) @@ -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. ! @@ -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) @@ -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) @@ -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 @@ -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 ) ] @@ -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) @@ -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)) @@ -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. @@ -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) diff --git a/components/eam/src/physics/cam/gw_drag.F90 b/components/eam/src/physics/cam/gw_drag.F90 index 1c6b7920e85a..99486f3c781b 100644 --- a/components/eam/src/physics/cam/gw_drag.F90 +++ b/components/eam/src/physics/cam/gw_drag.F90 @@ -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 !========================================================================== @@ -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 @@ -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 @@ -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. @@ -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", & @@ -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.